{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning how to compute time derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiation methods\n",
    "This notebook explores various choices of differentiation methods that can be used with PyKoopman's `KoopmanContinuous` class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optional parameter `differentiator` of the `KoopmanContinuous` class allows one to pass in custom time differentiation methods. `differentiator` should be callable with the call signature `differentiator(x, t)` where `x` is a 2D numpy `ndarray` with each example occupying a *row* and `t` is a 1D numpy `ndarray` containing the points in time for each row in `x`.\n",
    "\n",
    "Two common options for `differentiator` are 1) Methods from the [derivative](https://derivative.readthedocs.io/en/latest/) package, called via the `pykoopman.differentiation.Derivative` wrapper. 2) Entirely custom methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3.10.11'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('../src')\n",
    "\n",
    "from platform import python_version\n",
    "python_version() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-20T00:20:41.905651Z",
     "start_time": "2020-10-20T00:20:41.034185Z"
    }
   },
   "outputs": [],
   "source": [
    "from derivative import dxdt\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pykoopman as pk\n",
    "from pydmd import DMD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-20T00:20:42.078894Z",
     "start_time": "2020-10-20T00:20:41.907294Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHHCAYAAAB3K7g2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAx40lEQVR4nO3de3SU5YHH8d8kmBslA5iYi4a7iFwERYgByuUYReWwcnDbYMUgC9i6spaCF0JdEGkbtK6XKq2iIrSiYDWgW10UuehRIiCQlbjKgRgIahKMNhMIMULy7B80U8fcZsJc3nfm+zlnjubNM+88T94J88vzPheHMcYIAADAJqJCXQEAAABfEF4AAICtEF4AAICtEF4AAICtEF4AAICtEF4AAICtEF4AAICtEF4AAICtEF4AAICtEF4ABMX48eM1fvz4UFfDw3333SeHwxHqagDwEeEFgNvq1avlcDgUFxenL774otn3x48fr8GDB4egZgDwT4QXAM3U19dr+fLlfj3nW2+9pbfeesuv5zxb9957r+rq6kJdDQA+IrwAaGbYsGF6+umn9eWXX/rtnDExMYqJifHb+fyhU6dOiouLC3U1APiI8AKgmUWLFqmhocGr3pfTp09r2bJl6tu3r2JjY9WrVy8tWrRI9fX1HuVaGvPy+OOPa9CgQUpISFC3bt10+eWX64UXXpAkbdu2TQ6HQxs2bGj2mi+88IIcDocKCwtbrdepU6e0dOlSXXjhhYqLi9O5556rMWPGaPPmze4yLY15cTgcmjt3rjZu3KjBgwcrNjZWgwYN0qZNm5q9xhdffKF/+7d/U0pKirvcqlWr2v2ZATg7hBcAzfTu3Vu5uble9b7Mnj1bixcv1mWXXaZHHnlE48aNU35+vqZNm9bm855++mndcccdGjhwoB599FEtXbpUw4YN086dOyWdCTsZGRlau3Zts+euXbtWffv2VVZWVqvnv++++7R06VJNmDBBTzzxhH7961+rR48e2rt3b7vtf++99/Tv//7vmjZtmh588EF9++23uuGGG/T111+7y1RWVuqKK67Q22+/rblz5+qxxx5Tv379NGvWLD366KPtvgaAs2AA4B+ee+45I8ns3r3blJSUmE6dOpk77rjD/f1x48aZQYMGub8uKioykszs2bM9znPnnXcaSWbr1q0ezx03bpz76+uvv97jXC3Jy8szsbGxprq62n3s2LFjplOnTmbJkiVtPnfo0KFm0qRJbZZZsmSJ+eE/g5JMTEyMOXTokPvY//7v/xpJ5vHHH3cfmzVrlklLSzNVVVUez582bZpxOp3m5MmTbb42gI6j5wVAi/r06aObb75ZK1euVHl5eYtl3njjDUnS/PnzPY4vWLBAkvT666+3ev6uXbvq888/1+7du1stk5ubq/r6er388svuY+vXr9fp06c1ffr0NuvftWtXffzxxzp48GCb5VqSnZ2tvn37ur++5JJLlJiYqM8++0ySZIzRK6+8osmTJ8sYo6qqKvdj4sSJcrlcXvXwAOgYwguAVt177706ffp0q2Nfjhw5oqioKPXr18/jeGpqqrp27aojR460eu577rlHP/rRjzRy5EhdeOGFuv322/X+++97lBkwYIBGjBjhceto7dq1uuKKK5q95g/df//9qq6uVv/+/TVkyBDddddd+uijj9prsiSpR48ezY5169ZNf//73yVJX331laqrq7Vy5UolJyd7PGbOnClJOnbsmFevBcB3hBcArerTp4+mT5/eZu+LpA4t9HbxxRfrwIEDWrduncaMGaNXXnlFY8aM0ZIlSzzK5ebm6p133tHnn3+ukpISffDBB+32ukjS2LFjVVJSolWrVmnw4MF65plndNlll+mZZ55p97nR0dEtHjfGSJIaGxslSdOnT9fmzZtbfIwePbrd1wHQMZ1CXQEA1nbvvffq+eef1wMPPNDsez179lRjY6MOHjyoiy++2H28srJS1dXV6tmzZ5vn7ty5s3JycpSTk6PvvvtOU6dO1W9/+1vl5eW5pzBPmzZN8+fP14svvqi6ujqdc845ysnJ8aru3bt318yZMzVz5kydOHFCY8eO1X333afZs2f78BNoLjk5WV26dFFDQ4Oys7PP6lwAfEfPC4A29e3bV9OnT9dTTz2liooKj+9dd911ktRsds3DDz8sSZo0aVKr5/3+zB3pzDowAwcOlDFGp06dch9PSkrStddeq+eff15r167VNddco6SkpHbr/cPz/+hHP1K/fv2aTeHuiOjoaN1www165ZVXVFxc3Oz7X3311Vm/BoDW0fMCoF2//vWv9Ze//EUHDhzQoEGD3MeHDh2qGTNmaOXKlaqurta4ceO0a9curVmzRlOmTNGECRNaPefVV1+t1NRUjR49WikpKfrkk0/0xBNPaNKkSerSpYtH2dzcXP3rv/6rJGnZsmVe1XngwIEaP368hg8fru7du+vDDz/Uyy+/rLlz53bgJ9Dc8uXLtW3bNmVmZmrOnDkaOHCgvvnmG+3du1dvv/22vvnmG7+8DoDmCC8A2tWvXz9Nnz5da9asafa9Z555Rn369NHq1au1YcMGpaamKi8vr9nYlR/6+c9/rrVr1+rhhx/WiRMndMEFF+iOO+7Qvffe26zs5MmT1a1bNzU2Nupf/uVfvKrzHXfcoddee01vvfWW6uvr1bNnT/3mN7/RXXfd5V2j25GSkqJdu3bp/vvvV0FBgf74xz/q3HPP1aBBg1q8xQbAfxymaQQaAFjU6dOnlZ6ersmTJ+vZZ58NdXUAhBhjXgBY3saNG/XVV18pNzc31FUBYAH0vACwrJ07d+qjjz7SsmXLlJSUxMJvACTR8wLAwv70pz/ptttu03nnnac///nPoa4OAIug5wUAANgKPS8AAMBWCC8AAMBWwm6dl8bGRn355Zfq0qVLh/ZbAQAAwWeM0fHjx5Wenq6oqLb7VsIuvHz55ZfKyMgIdTUAAEAHHD16VBdccEGbZcIuvDQtK3706FElJiaGuDYAAMAbNTU1ysjIaLY9SEvCLrw03SpKTEwkvAAAYDPeDPlgwC4AALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgssq9xVpx0lVSp31YW6KgAACwm7jRkRHtbvLlNewX41GinKIeVPHaKcET1CXS0AgAXQ8wLLKXfVuYOLJDUaaVFBMT0wAABJhBdYUGlVrTu4NGkwRoerToamQgAASyG8wHJ6J3VWlMPzWLTDoV5JCaGpEADAUggvsJw0Z7zypw5RtONMgol2OPS7qYOV5owPcc0AAFbAgF1YUs6IHhrbP1mHq06qV1ICwQUA4EZ4gWWlOeMJLQCAZrhtBAAAbIXwAgAAbIXwAgAAbCWg4eXdd9/V5MmTlZ6eLofDoY0bN7b7nO3bt+uyyy5TbGys+vXrp9WrVweyirA5thAAgMgT0PBSW1uroUOHasWKFV6VLy0t1aRJkzRhwgQVFRVp3rx5mj17tt58881AVhM2tX53mUYv36qfPb1To5dv1frdZaGuEgAgCBzGGNN+MT+8kMOhDRs2aMqUKa2Wueeee/T666+ruLjYfWzatGmqrq7Wpk2bvHqdmpoaOZ1OuVwuJSYmnm21YVHlrjqNXr7VYyXeaIdD7y2cwAwlALAhXz6/LTXmpbCwUNnZ2R7HJk6cqMLCwhDVCFblyxYC3FoCgPBiqXVeKioqlJKS4nEsJSVFNTU1qqurU3x887+o6+vrVV9f7/66pqYm4PVE6DVtIfDDnpcfbiHgy+7U5a46lVbVqndSZ3pvAMDCLNXz0hH5+flyOp3uR0ZGRqirhCDwZgsBX3anZvwMANiHpXpeUlNTVVlZ6XGssrJSiYmJLfa6SFJeXp7mz5/v/rqmpoYAEyHa20KgrVtL3oScsf2T6YEBAAuyVHjJysrSG2+84XFs8+bNysrKavU5sbGxio2NDXTVYFFtbSHg7a0lb0MOAMAaAnrb6MSJEyoqKlJRUZGkM1Ohi4qKVFZ2pks+Ly9Pubm57vK/+MUv9Nlnn+nuu+/Wp59+qj/+8Y966aWX9Ktf/SqQ1USY8nZ36qaQ830thRwAgDUEtOflww8/1IQJE9xfN93emTFjhlavXq3y8nJ3kJGk3r176/XXX9evfvUrPfbYY7rgggv0zDPPaOLEiYGsJsKYN7tTN4WcRQXFajCm1ZADALCGoK3zEiys84KOKnfVtRlyAACB48vnt6XGvACh1Nb4GQCAddh+qjQAAIgshBcAAGArhBfAB/7aaoAtCwCg4xjzAnjJl60GgnEeAIhU9LwAXvBlq4FgnAcAIhnhBfCCL7tYB+M8ABDJCC+AF/y1Ci+r+QLA2SO8AF7wdquBYJ0HACIZK+wCPvDXKrys5gsAnlhhFwgQf63Cy2q+ANBx3DYCAAC2QniB37EAGwAgkLhtBL9iATYAQKDR8wK/YQE2AEAwEF7gNyzABgAIBsIL/IYF2AAAwUB4gd+wABsAIBgYsAu/yhnRQ2P7J7MAGwAgYAgvkHRmsG1pVa16J3VuNXB4U0bybgE2b88FAMAPEV7g1fRmf06BZjo1AOBsMOYlwnkzvdmfU6CZTg0AOFuElwjnzfRmf06BZjo1AOBsEV4inDfTm/05BZrp1ACAs0V4iXDeTG/25xRoplMDAM6Wwxhj2i9mHzU1NXI6nXK5XEpMTAx1dWyj3FXX7vRmb8r48/XsiplUAOA7Xz6/mW0ESd5Nb/amjD9fz46YSQUAgcdtI8BPmEkFAMFBeAH8hJlUABAchBfAT5hJBQDBQXgB/ISZVAAQHAzYBfyIjSkBIPAIL4CfhetMKgCwCm4bARZV7qrTjpIqZisBwA/Q8wJYkC/rxbAoHoBIQ3gBLKa19WLG9k9uFk5YFA9AJOK2EWAx3q4Xw6J4ACIV4QWwGG/Xi2FRPACRivACWIy368WwKB6ASMWYlwjAgE778Wa9mKaQs6igWA3GsCgegIhBeAlzDOi0L2/Wi2FRPACRiNtGYYwBnZEhzRmvrL7nElwARAzCSxhjQCcAIBwRXsIYAzoBAOGI8BLG2OUYTdhqAEA4YcBumGNAJxi0DSDc0PMSARjQGbkYtA0gHBFegDDGoG0A4YjwAoQxBm0DCEeEFyCMMWgbQDhiwC4Q5hi0DSDcEF6ACODNVgMAYBfcNgIAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALZCeAEAALYSlPCyYsUK9erVS3FxccrMzNSuXbtaLbt69Wo5HA6PR1xcXDCqCQAAbCDg4WX9+vWaP3++lixZor1792ro0KGaOHGijh071upzEhMTVV5e7n4cOXIk0NUEAAA2EfDw8vDDD2vOnDmaOXOmBg4cqCeffFIJCQlatWpVq89xOBxKTU11P1JSUgJdTSDilbvqtKOkSuWuulBXBQDaFNDw8t1332nPnj3Kzs7+5wtGRSk7O1uFhYWtPu/EiRPq2bOnMjIydP311+vjjz9utWx9fb1qamo8HgB8s353mUYv36qfPb1To5dv1frdZaGuEgC0KqDhpaqqSg0NDc16TlJSUlRRUdHicy666CKtWrVKr776qp5//nk1NjZq1KhR+vzzz1ssn5+fL6fT6X5kZGT4vR1AOCt31SmvYL8azZmvG420qKCYHhgAlmW52UZZWVnKzc3VsGHDNG7cOBUUFCg5OVlPPfVUi+Xz8vLkcrncj6NHjwa5xoC9lVbVuoNLkwZjdLjqZGgqBADt6BTIkyclJSk6OlqVlZUexysrK5WamurVOc455xxdeumlOnToUIvfj42NVWxs7FnXFYhUvZM6K8ohjwAT7XCoV1JC6CoFAG0IaM9LTEyMhg8fri1btriPNTY2asuWLcrKyvLqHA0NDdq/f7/S0tICVU0goqU545U/dYiiHQ5JZ4LL76YOVpozPsQ1A4CWBbTnRZLmz5+vGTNm6PLLL9fIkSP16KOPqra2VjNnzpQk5ebm6vzzz1d+fr4k6f7779cVV1yhfv36qbq6Wr///e915MgRzZ49O9BVtaVyV51Kq2rVO6kzHzbosJwRPTS2f7IOV51Ur6QE3ksALC3g4SUnJ0dfffWVFi9erIqKCg0bNkybNm1yD+ItKytTVNQ/O4D+/ve/a86cOaqoqFC3bt00fPhw7dixQwMHDgx0VW1n/e4y90DLKIeUP3WIckb0CHW1YFNpznhCCwBbcBhjTPvF7KOmpkZOp1Mul0uJiYmhrk7AlLvqNHr51mbjFN5bOIEPIACA7fjy+W252UbwDjNEAACRivBiU00zRL6PGSIAgEhAeLEpZogAACJVwAfsInCYIQIAiESEF5tjhggAINJw2wgAANgK4QUAANgK4QUAANgK4QUAANgK4QUAANgK4QUAANgK4QWA18pdddpRUqVyV12oqwIggrHOCwCvsIs5AKug5wVAu8pdde7gIkmNRlpUUEwPDICQILwAaBe7mAOwEsILgHaxizkAKyG8AGgXu5gDsBIG7ALwCruYA7AKwksIlLvqVFpVq95JnfkAgK2wizkAKyC8BBnTTQEAODuMeQkippsCAHD2CC9B5Mt0U1YyhV3x3gUQaNw2CqKm6abfDzAtTTfl1hLsivcugGCg5yWIvJluyq0l2JWv7116aAB0FD0vQdbedNO2bi0xywNW5st719seGmbmAWgJ4SUE2ppu6u2tJcBqvH3vttZDM7Z/ssfvBbegALSG20YWw0qmsCtv37veDFzn9imAttDzYkGsZAq78ua9600PDbdPAbSFnheLSnPGK6vvufxDDdtp773rTQ8NG0ECaAs9LwCCrr0emqaAs6igWA3GcPsUgAfCC4CQaG+fJG6fAmgN4QWAZbERJICWMOYFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuEFAADYCuHFz8pdddpRUqVyV12oqwIAQFhiV2k/Wr+7THkF+9VopCiHlD91iHJG9Ah1tQAACCv0vPhJuavOHVwkqdFIiwqK6YEBAMDPCC9+UlpV6w4uTRqM0eGqk6GpEAAAYYrw4ie9kzoryuF5LNrhUK+khNBUCACAMEV48ZM0Z7zypw5RtONMgol2OPS7qYOV5owPcc0AAAgvDNj1o5wRPTS2f7IOV51Ur6QEggsAAAFAePGzNGc8oQUAgADithEAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwgsAALAVwosP2DEaAIDQY50XL7FjNAAA1kDPixfYMRoAAOsgvHiBHaMBALAOwosX2DEaAADrILx4gR2jAQCwDgbseokdowEAsAbCiw/YMRoAgNDjthEAALAVwgsAALCVoISXFStWqFevXoqLi1NmZqZ27drVZvm//vWvGjBggOLi4jRkyBC98cYbwagmAACwgYCHl/Xr12v+/PlasmSJ9u7dq6FDh2rixIk6duxYi+V37NihG2+8UbNmzdK+ffs0ZcoUTZkyRcXFxYGuKgAAsAGHMca0X6zjMjMzNWLECD3xxBOSpMbGRmVkZOg//uM/tHDhwmblc3JyVFtbq7/97W/uY1dccYWGDRumJ598st3Xq6mpkdPplMvlUmJiov8aAgAAAsaXz++A9rx899132rNnj7Kzs//5glFRys7OVmFhYYvPKSws9CgvSRMnTmy1fH19vWpqajweAAAgfAU0vFRVVamhoUEpKSkex1NSUlRRUdHicyoqKnwqn5+fL6fT6X5kZGT4p/IAAMCSbD/bKC8vTy6Xy/04evRoqKsEAAACKKCL1CUlJSk6OlqVlZUexysrK5Wamtric1JTU30qHxsbq9jYWP9UGAAAWF5Ae15iYmI0fPhwbdmyxX2ssbFRW7ZsUVZWVovPycrK8igvSZs3b261PAAAiCwB3x5g/vz5mjFjhi6//HKNHDlSjz76qGprazVz5kxJUm5urs4//3zl5+dLkn75y19q3Lhx+q//+i9NmjRJ69at04cffqiVK1cGuqoAAMAGAh5ecnJy9NVXX2nx4sWqqKjQsGHDtGnTJveg3LKyMkVF/bMDaNSoUXrhhRd07733atGiRbrwwgu1ceNGDR48ONBVBQAANhDwdV6CjXVeAACwH8us8wIAAOBvhBcAAGArhBcAAGArhBcAAGArhBcAAGArhBcAAGArhBcAAOC1cleddpRUqdxVF7I6BHyROgAAEB7W7y5TXsF+NRopyiHlTx2inBE9gl4Pel4AAEC7yl117uAiSY1GWlRQHJIeGMILAFuzQhc2EAlKq2rdwaVJgzE6XHUy6HXhthEA27JKFzYQCXondVaUQx4BJtrhUK+khKDXhZ4XALZkpS5sIBKkOeOVP3WIoh0OSWeCy++mDlaaMz7odaHnBYAttdWFHYp/TIFIkDOih8b2T9bhqpPqlZQQst81wgsAW7JSFzYQSdKc8SH/A4HbRgBsyUpd2ACCi54XALZllS5sAMFFeAFga1bowgYQXNw2AgAAtkJ4AQAAkuyz6CO3jQAAgK0WfaTnBQCACGe3RR8JLwAARDgr7VvkDcILAAARrmnRx++z8qKPhBcAACKc3RZ9ZMAuAACw1aKPhBcAACDJPos+ctsIAADYCuEFQNizy8JbALzDbSMAYc1OC28B8A49LwDClt0W3gLgHcILgLBlt4W3AHiH8AIgbNlt4S0A3iG8AAhbdlt4C4B3GLALIKzZaeEtAN4hvAAIe3ZZeAuAd7htBAAAbIXwAgBAmAu3hRq5bQQAgI2Vu+pUWlWr3kmdW7w9Go4LNRJeAACwqfaCSWsLNY7tn2zrcWDcNgIAwIa8WUE6XBdqJLwAAGBD3gSTcF2okfACAIANeRNMwnWhRsa8AABgQ03BZFFBsRqMaTWYhONCjYQXAABsyttgEm4LNRJeAACwsXALJt5gzAsAALAVwgsAALAVwgsA+CDcllkH7IgxLwCg9pdYl8JzmXXAjggvACKeN6EkXJdZB+yI20YAIpo3S6xL4bvMOmBHhBcAEc3bUBKuy6wDdkR4ARDRvA0l4brMOmBHjHkBENG8XWJdCs9l1gE7IrwAiHi+hJJIXM0UsBrCCwCIUALYCWNeAACArRBeAACArRBeAACArRBeAACArRBeAACArRBeAACArRBeAACArRBeAACwqHJXnXaUVDXbKDTSsUgdAAAWtH53mXvH8yiHlD91iHJG9Ah1tSyBnhcAACym3FXnDi6S1GikRQXF9MD8A+EFAACLKa2qdQeXJg3G6HDVydBUyGIILwAAWEzvpM6Kcngei3Y41CspITQVspiAhpdvvvlGN910kxITE9W1a1fNmjVLJ06caPM548ePl8Ph8Hj84he/CGQ1AQCwlDRnvPKnDlG040yCiXY49Lupg9k89B8COmD3pptuUnl5uTZv3qxTp05p5syZuvXWW/XCCy+0+bw5c+bo/vvvd3+dkEDSBABElpwRPTS2f7IOV51Ur6QEgsv3BCy8fPLJJ9q0aZN2796tyy+/XJL0+OOP67rrrtNDDz2k9PT0Vp+bkJCg1NTUQFUNAABbSHPGE1paELDbRoWFheratas7uEhSdna2oqKitHPnzjafu3btWiUlJWnw4MHKy8vTyZOtD1Cqr69XTU2NxwMAAISvgPW8VFRU6LzzzvN8sU6d1L17d1VUVLT6vJ/97Gfq2bOn0tPT9dFHH+mee+7RgQMHVFBQ0GL5/Px8LV261K91BwAA1uVzeFm4cKEeeOCBNst88sknHa7Qrbfe6v7/IUOGKC0tTVdeeaVKSkrUt2/fZuXz8vI0f/5899c1NTXKyMjo8OsDAABr8zm8LFiwQLfcckubZfr06aPU1FQdO3bM4/jp06f1zTff+DSeJTMzU5J06NChFsNLbGysYmNjvT4fAACwN5/DS3JyspKTk9stl5WVperqau3Zs0fDhw+XJG3dulWNjY3uQOKNoqIiSVJaWpqvVQUAAGEoYAN2L774Yl1zzTWaM2eOdu3apffff19z587VtGnT3DONvvjiCw0YMEC7du2SJJWUlGjZsmXas2ePDh8+rNdee025ubkaO3asLrnkkkBVFQAA2EhAF6lbu3atBgwYoCuvvFLXXXedxowZo5UrV7q/f+rUKR04cMA9mygmJkZvv/22rr76ag0YMEALFizQDTfcoP/+7/8OZDUBAICNOIwxpv1i9lFTUyOn0ymXy6XExMRQVwcAAHjBl89v9jYCAAC2QngBAAC2QngBAAC2QngBAAC2QngBAAC2QngBAAC2QngBAD8rd9VpR0mVyl11oa4KEJYCtqs0AESi9bvLlFewX41GinJI+VOHKGdEj1BXCwgr9LwAgJ+Uu+rcwUWSGo20qKCYHhjAzwgvAOAnpVW17uDSpMEYHa46GZoKAWGK8AIAftI7qbOiHJ7Hoh0O9UpKCE2FgDBFeAEAP0lzxit/6hBFO84kmGiHQ7+bOlhpzvgQ1wwILwzYBQA/yhnRQ2P7J+tw1Un1SkoguKBV5a46lVbVqndSZ94nPiK8AICfpTnj+TBCm5iVdna4bQQAQBAxK+3sEV4AAAgiZqWdPcILAABBxKy0s0d4AQAgiJiVdvYYsAsAQJAxK+3sEF4AAPAzb6ZBMyut4wgvAADIf+uuMA068AgvAICI56/A0do06LH9k+ll8SMG7AJACJS76rSjpIq1PSzAn+uuMA06OOh5AYAg8/avfJaPD462AoevP/emadDfPx/ToP2PnhcACCJv/8pfv7tMo5dv1c+e3qnRy7dq/e6yENQ2Mvhz3RWmQQcHPS8AEETe/JXPuIngagociwqK1WDMWQcOpkEHHuEFAILIm9sK/ryNAe/4O3AwDTqwuG0EAEHkzW0Flo8PjTRnvLL6nkvosAF6XgAgyNr7K9/ftzGAcEN4AYAQaO+2AuMmgNYRXgDAohg3AbSMMS8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AAMBWCC8AEAHKXXXaUVKlclddqKsCnDW2BwCAMLd+d5nyCvar0UhRDil/6hDljOgR6moBHUbPCwCEsXJXnTu4SFKjkRYVFNMDA1sjvABAGCutqnUHlyYNxuhw1cnQVAjwA8ILAISx3kmdFeXwPBbtcKhXUkJoKgT4AeEFAMJYmjNe+VOHKNpxJsFEOxz63dTBSnPGh7hmQMcxYBcAwlzOiB4a2z9Zh6tOqldSQkQGl3JXnUqratU7qfNZtd9f58HZIbwAQARIc8ZH7Ietv2ZbMWvLOrhtBAAIW/6abcWsLWshvAAAwpa/Zlsxa8taCC8AgLDlr9lWzNqyFsILAMDW2tr6wF+zrZi1ZS0OY4xpv5h91NTUyOl0yuVyKTExMdTVAQAEkLeDaMtddX6ZbeWv86A5Xz6/mW0EALCl1gbRju2f3CxY+Gu2VSTP2rISbhsBAGyJQbSRi/ACALAlBtFGLsILAMCWGEQbuRjzAgCwrPaW42frg8hEeAEASLLevj3eziRiEG3kIbwAACy3b48vM4kQeRjzAgARzor79jCTCG0hvABAhLNiUGAmEdpCeAGACBeKoNDWkv4SM4nQNsa8AECEawoKiwqK1WBMwIOCt+NrmEmE1rC3EQDYmD9nCAVj355yV51GL9/qcZsq2uHQewsnEE4iHHsbAUAE8PcMIW+mHJ9tWGprfA3hBd4ivACADYViKrE/wlLT+Jof9rwwEBe+YMAuANhQsGcI+Ws6NQNx4Q8BCy+//e1vNWrUKCUkJKhr165ePccYo8WLFystLU3x8fHKzs7WwYMHA1VFALCtYM8Q8mdYyhnRQ+8tnKAX51yh9xZOCOlieLCngIWX7777Tj/5yU902223ef2cBx98UH/4wx/05JNPaufOnercubMmTpyob7/9NlDVBABbCnYPhr/DUpozXll9z6XHBR0S8NlGq1ev1rx581RdXd1mOWOM0tPTtWDBAt15552SJJfLpZSUFK1evVrTpk3z6vWYbQQgkgRjhlCT9bvLmk2nptcE/mLL2UalpaWqqKhQdna2+5jT6VRmZqYKCwtbDS/19fWqr693f11TUxPwugKAVQRzU0LWXYFVWGbAbkVFhSQpJSXF43hKSor7ey3Jz8+X0+l0PzIyMgJaTwCIZNzugRX4FF4WLlwoh8PR5uPTTz8NVF1blJeXJ5fL5X4cPXo0qK8PAACCy6fbRgsWLNAtt9zSZpk+ffp0qCKpqamSpMrKSqWlpbmPV1ZWatiwYa0+LzY2VrGxsR16TQAAYD8+hZfk5GQlJycHpCK9e/dWamqqtmzZ4g4rNTU12rlzp08zlgAAQHgL2JiXsrIyFRUVqaysTA0NDSoqKlJRUZFOnDjhLjNgwABt2LBBkuRwODRv3jz95je/0Wuvvab9+/crNzdX6enpmjJlSqCqCQAAbCZgs40WL16sNWvWuL++9NJLJUnbtm3T+PHjJUkHDhyQy+Vyl7n77rtVW1urW2+9VdXV1RozZow2bdqkuLi4QFUTAOBn/twsEmgJu0oDAPzG35tFInL48vltmanSAAB789f+R0B7CC8AAL8I9maRiFyEFwCAXwR7s0hELsILAMAvgr1ZJCKXZfY2AgDYH/sfIRgILwAAvwrmZpGITNw2AgAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAthJ2exsZYyRJNTU1Ia4JAADwVtPndtPneFvCLrwcP35ckpSRkRHimgAAAF8dP35cTqezzTIO403EsZHGxkZ9+eWX6tKlixwOh1/PXVNTo4yMDB09elSJiYl+PbcVhHv7pPBvI+2zv3BvI+2zv0C10Rij48ePKz09XVFRbY9qCbuel6ioKF1wwQUBfY3ExMSwfVNK4d8+KfzbSPvsL9zbSPvsLxBtbK/HpQkDdgEAgK0QXgAAgK0QXnwQGxurJUuWKDY2NtRVCYhwb58U/m2kffYX7m2kffZnhTaG3YBdAAAQ3uh5AQAAtkJ4AQAAtkJ4AQAAtkJ4AQAAthLR4WXFihXq1auX4uLilJmZqV27drVZ/q9//asGDBiguLg4DRkyRG+88YbH940xWrx4sdLS0hQfH6/s7GwdPHgwkE1oly9tfPrpp/XjH/9Y3bp1U7du3ZSdnd2s/C233CKHw+HxuOaaawLdjFb50r7Vq1c3q3tcXJxHGbtfw/Hjxzdro8Ph0KRJk9xlrHQN3333XU2ePFnp6elyOBzauHFju8/Zvn27LrvsMsXGxqpfv35avXp1szK+/m4Hiq/tKygo0FVXXaXk5GQlJiYqKytLb775pkeZ++67r9n1GzBgQABb0Tpf27d9+/YW358VFRUe5axy/STf29jS75fD4dCgQYPcZaxyDfPz8zVixAh16dJF5513nqZMmaIDBw60+zwrfBZGbHhZv3695s+fryVLlmjv3r0aOnSoJk6cqGPHjrVYfseOHbrxxhs1a9Ys7du3T1OmTNGUKVNUXFzsLvPggw/qD3/4g5588knt3LlTnTt31sSJE/Xtt98Gq1kefG3j9u3bdeONN2rbtm0qLCxURkaGrr76an3xxRce5a655hqVl5e7Hy+++GIwmtOMr+2TzqwI+f26HzlyxOP7dr+GBQUFHu0rLi5WdHS0fvKTn3iUs8o1rK2t1dChQ7VixQqvypeWlmrSpEmaMGGCioqKNG/ePM2ePdvjA74j74tA8bV97777rq666iq98cYb2rNnjyZMmKDJkydr3759HuUGDRrkcf3ee++9QFS/Xb62r8mBAwc86n/eeee5v2el6yf53sbHHnvMo21Hjx5V9+7dm/0OWuEavvPOO7r99tv1wQcfaPPmzTp16pSuvvpq1dbWtvocy3wWmgg1cuRIc/vtt7u/bmhoMOnp6SY/P7/F8j/96U/NpEmTPI5lZmaan//858YYYxobG01qaqr5/e9/7/5+dXW1iY2NNS+++GIAWtA+X9v4Q6dPnzZdunQxa9ascR+bMWOGuf766/1d1Q7xtX3PPfeccTqdrZ4vHK/hI488Yrp06WJOnDjhPmala/h9ksyGDRvaLHP33XebQYMGeRzLyckxEydOdH99tj+zQPGmfS0ZOHCgWbp0qfvrJUuWmKFDh/qvYn7iTfu2bdtmJJm///3vrZax6vUzpmPXcMOGDcbhcJjDhw+7j1n1Gh47dsxIMu+8806rZazyWRiRPS/fffed9uzZo+zsbPexqKgoZWdnq7CwsMXnFBYWepSXpIkTJ7rLl5aWqqKiwqOM0+lUZmZmq+cMpI608YdOnjypU6dOqXv37h7Ht2/frvPOO08XXXSRbrvtNn399dd+rbs3Otq+EydOqGfPnsrIyND111+vjz/+2P29cLyGzz77rKZNm6bOnTt7HLfCNeyI9n4P/fEzs5LGxkYdP3682e/gwYMHlZ6erj59+uimm25SWVlZiGrYMcOGDVNaWpquuuoqvf/+++7j4Xb9pDO/g9nZ2erZs6fHcSteQ5fLJUnN3m/fZ5XPwogML1VVVWpoaFBKSorH8ZSUlGb3XptUVFS0Wb7pv76cM5A60sYfuueee5Senu7xJrzmmmv05z//WVu2bNEDDzygd955R9dee60aGhr8Wv/2dKR9F110kVatWqVXX31Vzz//vBobGzVq1Ch9/vnnksLvGu7atUvFxcWaPXu2x3GrXMOOaO33sKamRnV1dX5531vJQw89pBMnTuinP/2p+1hmZqZWr16tTZs26U9/+pNKS0v14x//WMePHw9hTb2TlpamJ598Uq+88opeeeUVZWRkaPz48dq7d68k//y7ZSVffvml/ud//qfZ76AVr2FjY6PmzZun0aNHa/Dgwa2Ws8pnYdjtKg3/WL58udatW6ft27d7DGqdNm2a+/+HDBmiSy65RH379tX27dt15ZVXhqKqXsvKylJWVpb761GjRuniiy/WU089pWXLloWwZoHx7LPPasiQIRo5cqTHcTtfw0jywgsvaOnSpXr11Vc9xoRce+217v+/5JJLlJmZqZ49e+qll17SrFmzQlFVr1100UW66KKL3F+PGjVKJSUleuSRR/SXv/wlhDULjDVr1qhr166aMmWKx3ErXsPbb79dxcXFIRs/5auI7HlJSkpSdHS0KisrPY5XVlYqNTW1xeekpqa2Wb7pv76cM5A60sYmDz30kJYvX6633npLl1xySZtl+/Tpo6SkJB06dOis6+yLs2lfk3POOUeXXnqpu+7hdA1ra2u1bt06r/4hDNU17IjWfg8TExMVHx/vl/eFFaxbt06zZ8/WSy+91KyL/oe6du2q/v372+L6tWTkyJHuuofL9ZPOzLhZtWqVbr75ZsXExLRZNtTXcO7cufrb3/6mbdu26YILLmizrFU+CyMyvMTExGj48OHasmWL+1hjY6O2bNni8Zf592VlZXmUl6TNmze7y/fu3VupqakeZWpqarRz585WzxlIHWmjdGaU+LJly7Rp0yZdfvnl7b7O559/rq+//lppaWl+qbe3Otq+72toaND+/fvddQ+XayidmcpYX1+v6dOnt/s6obqGHdHe76E/3heh9uKLL2rmzJl68cUXPaa4t+bEiRMqKSmxxfVrSVFRkbvu4XD9mrzzzjs6dOiQV39AhOoaGmM0d+5cbdiwQVu3blXv3r3bfY5lPgv9NvTXZtatW2diY2PN6tWrzf/93/+ZW2+91XTt2tVUVFQYY4y5+eabzcKFC93l33//fdOpUyfz0EMPmU8++cQsWbLEnHPOOWb//v3uMsuXLzddu3Y1r776qvnoo4/M9ddfb3r37m3q6uqC3j5jfG/j8uXLTUxMjHn55ZdNeXm5+3H8+HFjjDHHjx83d955pyksLDSlpaXm7bffNpdddpm58MILzbfffmv59i1dutS8+eabpqSkxOzZs8dMmzbNxMXFmY8//thdxu7XsMmYMWNMTk5Os+NWu4bHjx83+/btM/v27TOSzMMPP2z27dtnjhw5YowxZuHChebmm292l//ss89MQkKCueuuu8wnn3xiVqxYYaKjo82mTZvcZdr7mVm5fWvXrjWdOnUyK1as8PgdrK6udpdZsGCB2b59uyktLTXvv/++yc7ONklJSebYsWOWb98jjzxiNm7caA4ePGj2799vfvnLX5qoqCjz9ttvu8tY6foZ43sbm0yfPt1kZma2eE6rXMPbbrvNOJ1Os337do/328mTJ91lrPpZGLHhxRhjHn/8cdOjRw8TExNjRo4caT744AP398aNG2dmzJjhUf6ll14y/fv3NzExMWbQoEHm9ddf9/h+Y2Oj+c///E+TkpJiYmNjzZVXXmkOHDgQjKa0ypc29uzZ00hq9liyZIkxxpiTJ0+aq6++2iQnJ5tzzjnH9OzZ08yZMydk/6gY41v75s2b5y6bkpJirrvuOrN3716P89n9GhpjzKeffmokmbfeeqvZuax2DZumzv7w0dSmGTNmmHHjxjV7zrBhw0xMTIzp06ePee6555qdt62fWTD52r5x48a1Wd6YM1PD09LSTExMjDn//PNNTk6OOXToUHAb9g++tu+BBx4wffv2NXFxcaZ79+5m/PjxZuvWrc3Oa5XrZ0zH3qPV1dUmPj7erFy5ssVzWuUattQuSR6/U1b9LHT8owEAAAC2EJFjXgAAgH0RXgAAgK0QXgAAgK0QXgAAgK0QXgAAgK0QXgAAgK0QXgAAgK0QXgAAgK0QXgDYxvjx4zVv3rxQVwNAiBFeAACArbA9AABbuOWWW7RmzRqPY6WlperVq1doKgQgZAgvAGzB5XLp2muv1eDBg3X//fdLkpKTkxUdHR3imgEItk6hrgAAeMPpdComJkYJCQlKTU0NdXUAhBBjXgAAgK0QXgAAgK0QXgDYRkxMjBoaGkJdDQAhRngBYBu9evXSzp07dfjwYVVVVamxsTHUVQIQAoQXALZx5513Kjo6WgMHDlRycrLKyspCXSUAIcBUaQAAYCv0vAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFshvAAAAFv5f+4RLl7PaT30AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.linspace(0, 2, 50)\n",
    "x = np.sin(np.pi * t) + 0.1 * np.random.standard_normal(t.shape)\n",
    "x_dot = np.pi * np.cos(np.pi *  t)\n",
    "\n",
    "plt.plot(t, x, '.')\n",
    "plt.xlabel('t')\n",
    "plt.title('Noisy sine');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivative package\n",
    "\n",
    "All of the robust differentiation methods in the [derivative](https://derivative.readthedocs.io/en/latest/) package are available with the `pykoopman.differentiation.Derivative` wrapper class. One need only pass in the same keyword arguments to `Derivative` that one would pass to [derivative.dxdt](https://derivative.readthedocs.io/en/latest/api.html#dxdt-functional-interface).\n",
    "\n",
    "For example, we'll compute the derivative of a noisy sine function with `dxdt` and `Derivative` using spline-based numerical differentiation and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-20T00:20:42.281295Z",
     "start_time": "2020-10-20T00:20:42.081151Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAHHCAYAAACC1TOGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAACIsElEQVR4nOzdd3hUVfrA8e+dmj7pjYTQSQi9g3SQKoodVAR73bX3VXTVtev6c+0FXBWxIOqK9KoISO+ElpAACek9mUy5vz8iA5ckQCDJTJL38zzz6Jx77r3vzQwz75x7iqKqqooQQgghhBvo3B2AEEIIIZovSUSEEEII4TaSiAghhBDCbSQREUIIIYTbSCIihBBCCLeRREQIIYQQbiOJiBBCCCHcRhIRIYQQQriNJCJCCCGEcBtJRIRoYIqi8Oyzz9bb8WfNmoWiKKSkpNTbOTzN/v37GT16NBaLBUVR+PHHH90dUo1WrlyJoih8//339X6u6dOn06pVq3o/jxAXQhIR0eyc+KKu6bFu3Tp3hyhqadq0aezYsYMXX3yRL774gt69e7s7JGbPns2///1vd4chhMczuDsAIdzln//8J61bt65S3q5du3o9b1lZGQaD/NOrK2VlZaxdu5annnqKe++9193huMyePZudO3dy//33uzsUITyafBqKZmvcuHFu+eXs5eXV4OdsyrKysgAIDAx0byBCiPMit2aEOIP8/HymT5+OxWIhMDCQadOmsXXrVhRFYdasWa56w4YNY9iwYVX2r+4e/al9RL7//nsURWHVqlVV9v3www9RFIWdO3cCsH37dqZPn06bNm3w8vIiMjKSm2++mZycnHO6lgULFjB48GB8fX3x9/dnwoQJ7Nq1q0q8fn5+HD16lEmTJuHn50dYWBgPP/wwDodDU9fpdPL222/TpUsXvLy8CAsLY+zYsWzcuFFT78svv6RXr154e3sTHBzM5MmTSUtLO6eYt2zZwrhx4wgICMDPz4+RI0dqbp09++yzxMXFAfDII4+gKMoZ+0Sc6J/x7bff8txzz9GiRQv8/f256qqrKCgowGq1cv/99xMeHo6fnx833XQTVqu1ynHOdk3Dhg1j/vz5HD582HXL7/S4nE4nL774IjExMXh5eTFy5EgOHDhQ5Vzfffed61yhoaHccMMNHD16tEq9H3/8kc6dO+Pl5UXnzp2ZN29etX+DOXPm0KtXL/z9/QkICKBLly68/fbbNf7NhKhv0iIimq2CggKys7M1ZYqiEBISAoCqqlx22WX8/vvv3HnnnSQkJDBv3jymTZtWZzFMmDABPz8/vv32W4YOHarZ9s0335CYmEjnzp0BWLJkCYcOHeKmm24iMjKSXbt28dFHH7Fr1y7WrVuHoig1nueLL75g2rRpjBkzhldeeYXS0lLef/99Bg0axJYtWzRfkg6HgzFjxtCvXz9ef/11li5dyhtvvEHbtm256667XPVuueUWZs2axbhx47j11lux2+389ttvrFu3ztXS9OKLL/L0009zzTXXcOutt5KVlcU777zDkCFD2LJlyxlbMXbt2sXgwYMJCAjg0UcfxWg08uGHHzJs2DBWrVpFv379uOKKKwgMDOSBBx5gypQpjB8/Hj8/v7P+3V966SW8vb15/PHHOXDgAO+88w5GoxGdTkdeXh7PPvss69atY9asWbRu3ZpnnnnGte+5XNNTTz1FQUEBR44c4a233gKoEtfLL7+MTqfj4YcfpqCggFdffZXrr7+e9evXu+rMmjWLm266iT59+vDSSy9x/Phx3n77bdasWaP5+y1evJgrr7ySTp068dJLL5GTk8NNN91ETEyM5pxLlixhypQpjBw5kldeeQWAPXv2sGbNGu67776z/t2EqBeqEM3MzJkzVaDah9lsdtX78ccfVUB99dVXXWV2u10dPHiwCqgzZ850lQ8dOlQdOnRolXNNmzZNjYuL05QB6owZM1zPp0yZooaHh6t2u91Vlp6erup0OvWf//ynq6y0tLTK8b/++msVUFevXl3l+pKTk1VVVdWioiI1MDBQve222zT7ZmRkqBaLRVM+bdo0FdCcV1VVtUePHmqvXr1cz5cvX64C6t///vcqMTmdTlVVVTUlJUXV6/Xqiy++qNm+Y8cO1WAwVCk/3aRJk1STyaQePHjQVXbs2DHV399fHTJkiKssOTlZBdTXXnvtjMdTVVVdsWKFCqidO3dWKyoqXOVTpkxRFUVRx40bp6k/YMAAzetXm2uaMGFCldf+1BgSEhJUq9XqKn/77bdVQN2xY4eqqqpaUVGhhoeHq507d1bLyspc9X755RcVUJ955hlXWffu3dWoqCg1Pz/fVbZ48WIV0MRw3333qQEBAZr3mhDuJrdmRLP17rvvsmTJEs1jwYIFru2//vorBoNB0wqg1+v529/+VqdxXHvttWRmZrJy5UpX2ffff4/T6eTaa691lXl7e7v+v7y8nOzsbPr37w/A5s2bazz+kiVLyM/PZ8qUKWRnZ7seer2efv36sWLFiir73HnnnZrngwcP5tChQ67nc+fORVEUZsyYUWXfEy0zP/zwA06nk2uuuUZz3sjISNq3b1/teU9wOBwsXryYSZMm0aZNG1d5VFQU1113Hb///juFhYU17n82N954I0aj0fW8X79+qKrKzTffrKnXr18/0tLSsNvtF3xNp7vpppswmUyu54MHDwZw/Z03btxIZmYmd999t6Zf0YQJE4iPj2f+/PkApKens3XrVqZNm4bFYnHVu/jii+nUqZPmnIGBgZSUlLBkyZJzjlOI+ia3ZkSz1bdv3zN2Vj18+DBRUVFVmtQ7duxYp3GMHTsWi8XCN998w8iRI4HK2zLdu3enQ4cOrnq5ubk899xzzJkzh8zMTM0xCgoKajz+/v37ARgxYkS12wMCAjTPT/T3OFVQUBB5eXmu5wcPHiQ6Oprg4OAznldVVdq3b1/t9lMTgdNlZWVRWlpa7d86ISEBp9NJWloaiYmJNR7jTFq2bKl5fuILPDY2tkq50+mkoKCAkJCQC7qms8UQFBQE4Po7Hz58GKj+/RYfH8/vv/+uqVddTB07dtQkqXfffTfffvst48aNo0WLFowePZprrrmGsWPHnnPcQtQ1SUSEqAOKoqCqapXy0zt4VsdsNjNp0iTmzZvHe++9x/Hjx1mzZg3/+te/NPWuueYa/vjjDx555BG6d++On58fTqeTsWPH4nQ6azz+iW1ffPEFkZGRVbafPpRYr9efNeZz4XQ6URSFBQsWVHvMc+nLUV9qusaayk+8tnV5TWc7V30IDw9n69atLFq0iAULFrBgwQJmzpzJjTfeyOeff15v5xXiTCQREaIGcXFxLFu2jOLiYs0XTFJSUpW6QUFBmlsXJ5z4tXo21157LZ9//jnLli1jz549qKqquS2Tl5fHsmXLeO655zQdJ0+0dpxJ27ZtgcovoVGjRp1TPOdyzEWLFpGbm1tjq0jbtm1RVZXWrVtrWnbORVhYGD4+PtX+rffu3YtOp6vSetEQanNNZ+o8fC5OjAZKSkqq0pqVlJTk2n7iv9W9F6r7+5lMJiZOnMjEiRNxOp3cfffdfPjhhzz99NP1PoeOENWRPiJC1GD8+PHY7Xbef/99V5nD4eCdd96pUrdt27bs3bvXNacFwLZt21izZs05nWvUqFEEBwfzzTff8M0339C3b1/NZGsnfj2f/mv5XGbuHDNmDAEBAfzrX//CZrNV2X5qzOfqyiuvRFVVnnvuuSrbTsR4xRVXoNfree6556rErarqGYcd6/V6Ro8ezU8//aSZqv748ePMnj2bQYMGVbml1BBqc02+vr5nvGV2Nr179yY8PJwPPvhAM4R4wYIF7NmzhwkTJgCV/Wa6d+/O559/rjnfkiVL2L17t+aYp//NdTodXbt2Bah2mLIQDUFaRESztWDBAvbu3VulfODAgbRp04aJEydy0UUX8fjjj5OSkkKnTp344Ycfqv1yufnmm3nzzTcZM2YMt9xyC5mZmXzwwQckJiaeU6dKo9HIFVdcwZw5cygpKeH111/XbA8ICGDIkCG8+uqr2Gw2WrRoweLFi0lOTj7rsQMCAnj//feZOnUqPXv2ZPLkyYSFhZGamsr8+fO56KKL+M9//nPW45xq+PDhTJ06lf/7v/9j//79rttDv/32G8OHD+fee++lbdu2vPDCCzzxxBOkpKQwadIk/P39SU5OZt68edx+++08/PDDNZ7jhRdeYMmSJQwaNIi7774bg8HAhx9+iNVq5dVXX61VvHWlNtfUq1cvvvnmGx588EH69OmDn58fEydOPOdzGY1GXnnlFW666SaGDh3KlClTXMN3W7VqxQMPPOCq+9JLLzFhwgQGDRrEzTffTG5uLu+88w6JiYkUFxe76t16663k5uYyYsQIYmJiOHz4MO+88w7du3cnISGh7v5QQtRGww/UEcK9zjR8l9OG5ebk5KhTp05VAwICVIvFok6dOlXdsmVLlXqqqqpffvml2qZNG9VkMqndu3dXFy1adE7Dd09YsmSJCqiKoqhpaWlVth85ckS9/PLL1cDAQNVisahXX321euzYsSrHO3347gkrVqxQx4wZo1osFtXLy0tt27atOn36dHXjxo2uOtOmTVN9fX2rnHvGjBnq6R8Xdrtdfe2119T4+HjVZDKpYWFh6rhx49RNmzZp6s2dO1cdNGiQ6uvrq/r6+qrx8fHqPffcoyYlJVU5z+k2b96sjhkzRvXz81N9fHzU4cOHq3/88YemzvkM3/3uu+805Sf+Zhs2bKj2urOysmp9TcXFxep1112nBgYGaobR1hTDies4/X31zTffqD169FDNZrMaHBysXn/99eqRI0eqXNvcuXPVhIQE1Ww2q506dVJ/+OGHKu+/77//Xh09erQaHh6umkwmtWXLluodd9yhpqenn/VvJ0R9UVS1HntGCdEEpaSk0Lp1a2bOnMn06dPdHY4QQjRq0kdECCGEEG4jiYgQQggh3EYSESGEEEK4jfQREUIIIYTbSIuIEEIIIdxGEhEhhBBCuI1HT2jmdDo5duwY/v7+FzxdshBCCCEahqqqFBUVER0djU535jYPj05Ejh075pb1JIQQQghx4dLS0oiJiTljHY9ORPz9/YHKC3HHuhJCCCGEqL3CwkJiY2Nd3+Nn4tGJyInbMQEBAZKICCGEEI3MuXSrkM6qQgghhHAbSUSEEEII4TaSiAghhBDCbSQREUIIIYTbSCIihBBCCLeRREQIIYQQbiOJiBBCCCHcRhIRIYQQQriNJCJCCCGEcBtJRIQQQgjhNpKICCGEEMJtPHqtmfrgdKqsPZRDRICZSIs3fuZm9ycQQgghsDucZBZZOV5YTrCvibgQX7fE0ey+hXNKKvjs92TXc1+zwZWURASYiQzwIuKvh8kgDUZCCCEaL1VVySu1kVFQzvHCykd6QTmZReVkFVlR1cp6IxMiJBFpKMcLyzXPS6x2DmXZOZRVoilXFAj2NREb5MOQDmF0jbGc0yqCQgghhDvlFFtZvjeTXccKOV5YToXdedZ9Mk77bmxIzS4RySg4tz+2qkJOcQU5xRVsTcsnwuLF6E4RDGwbKi0lQgghPE5ydgmLd2WwISUP9URTxznKKCirp6jOrtklIga9QlSgF5mFVhzOc3+hjheU88Xaw8zbcpThHcMZHh+OxdtYj5EKIYQQZ+Z0qmw7ks+iXcfZf7zovI5h8TYS5m9GVVW3tPwram3TpgZUWFiIxWKhoKCAgICAOj22w6mSU2zleKGV9IIyjhdZySgo43ihlbySirPur9cpDGgbwujESFoEetdpbEIIIcSZlNscrD2Yw+LdGWQWWs9a32zUEe7vRZSlsg9kpOVEf0gzPqa6b5Oozfd3vSYi77//Pu+//z4pKSkAJCYm8swzzzBu3Lhz2r8+E5EzKbc5yCy0cii7mBV7MzmSd+Ymq84tLFzeowWtQt3T0UcIIUTzYLU7WLAjg+V7Mymx2musZzLoGNQ+lJ4tg4gM8CLQx9igrR0ek4j873//Q6/X0759e1RV5fPPP+e1115jy5YtJCYmnnV/dyUip1JVld3phSzadZxdRwtqrKfTKVzZM4YxiRHSqVUIIUSdO5pfxgcrD3Isv+YfxxYfI6MSIhjaIQxfN05P4TGJSHWCg4N57bXXuOWWW85a1xMSkVMdyStlye7jrD2YU2P/ks4tLNwyuDUBXtJ/RAghxIVTVZXf9mcze30qNkf1I2Big30YnRhB31bBGPTuH1BRm+/vBkuXHA4H3333HSUlJQwYMKChTlunYoJ8uOmi1lzRM4YVezOrbRrbebSAZ3/exe1D2hAf6f7kSQghRONVVuHgv2tT+DM5t9rtXWMCGZ0YQXykf6Ntja/3FpEdO3YwYMAAysvL8fPzY/bs2YwfP77aularFav1ZKebwsJCYmNjPaZF5HQVdie/bD/GrzvSOf2vqChwSddoLu0WjU7XON8cQggh3Cc5u4SPVh+stjNqbLAPtwxqTWywjxsiOzuPujVTUVFBamoqBQUFfP/993zyySesWrWKTp06Van77LPP8txzz1Up99RE5IRdxwr45LdkCstsVbZ1iPTn9sFtCPI1uSEyIYQQjY2qqizZfZzvNx2pthvA8Phwrukd69FzWnlUInK6UaNG0bZtWz788MMq2xpbi8ipCspsfPrbIXYdK6yyzdds4JZBrekWG9jwgQkhhGg0iq12Pv0tme1H8qts8zbpuemiVvSKC274wGrJI/uInOB0OjXJxqnMZjNms7mBI6obFm8jD1zcgQU7M/hh81HNrHYlVjv/t2w/YzpHcnWvmEZ7H08IIUT9OZRVzHsrD1Y7l1WbMF/uGNqWUL/G+R15JvWaiDzxxBOMGzeOli1bUlRUxOzZs1m5ciWLFi2qz9O6jaIojO8SRYcIPz5cdYjc095Mi3ZmUG5zMLV/nCQjQgghXA5mFfPm4n2U2xxVto3rEsWk7tEeMRqmPtRrIpKZmcmNN95Ieno6FouFrl27smjRIi6++OL6PK3btQv3Z8alicxak8yW1HzNtlVJWegUhev7tZRkRAghBMnZJby1pGoS4u9l4NbBbejcwuKmyBpGs53ivSGoqsryvZnM2ZCG87QOR6M6RTC5T6wkI0II0YwdzinhtUVJlFVok5AOkf7cMaQNgT6Nc6BDbb6/m2Y7j4dQFIWRCRHcObRNlYRj6e7jfLsxrdYrJAohhGga0nJLeWPxvipJSOcWFh4Y1aHRJiG1JYlIA+gVF8wd1SQji3dVDs+SZEQIIZqXtNxSXluUVGVSzE7RAdwzvJ1HD82ta83nSt2sT6tgbhvcmtPvxCzcmcG8LUclGRFCiGbiWH4ZbyyumoTER/lz74jmlYSAJCINql+bEG4eVDUZmb89nZ+3HXNPUEIIIRpMekEZry1Koqhcm4R0iPTnbyPaYzbo3RSZ+0gi0sAGtg3lpouqJiM/bz0myYgQQjRhGQXlvLYoqcos3O3C/bhvZHu8jM0vCQFJRNzionah3DigVZXyn7YcZf729IYPSAghRL3KLCzn1UV7KSjVJiFtw/144OIOzTYJAUlE3GZIhzCmDoirUv7D5iOs3pflhoiEEELUh2KrndcXJ1VJQlqH+nL/qObbEnKCJCJuNKxjONf3b1ml/Mt1h0nOLnFDREIIIeqS06ny4aqD5BRrZ9qOC/HlwdEd8DE1+EorHkcSETcbER/BlL7aZMThVHl3xQGKyquu5iuEEKLx+HHrUXafthhqbLCPJCGnkETEA4zqFMHFnSI0ZXklFXy0+lCVGVmFEEI0DltS86r0+7N4G7l/VHv8zJKEnCCJiIe4qlcM7SP8NWW7jxUyb8tRN0UkhBDifB0vLOeT35M1ZTqdwl3D2jabGVPPVbNORL5b/ADL/3iVpP3zKS5y79BZg17HXUPbYvExasp/3ZHO5tQ8N0UlhBCitsptDt5dcYDy06Zuv6Z3bJUfnPVy/rI8du/9kWV/vMIff7wGNUyYeSh5Gcv/eLXe4zmbZts2VFR4lFUZ6yqfHJgHwNM+HYkIbgeBcRDYsvJhiQG98QxHqjsWHyN3D2vLKwuTNLdkPv09megJ3kRavBokDiGEEOdHVVX+uzaFo3llmvK+rYMZlRBePyetKIGsvXB8Nz8eXsyyooOoVH6HRCtGBva7v9rvseLSbA7nJtVPTLXQbBORY8e3ap4bUQizlkLGjsrHCXojxPbjaGQnoloNQ6er32FW7cL9ubZ3LF//meoqK6+ozK6fmpDQ7Id5CSGEJ1u2J5P1h3I1ZS2CvJk2sFXdrbZuLYbM3X899kDeYfgr8bA4i1xJCIAdwGGrNhEx6M3Yne4fFNFsE5GjWbs0zyMVI7rq3iQOG6XJq3ntwNf4rfOhX1Q/+iVOITy8c73FNjIhnOTsEtYdynGVHcsvY9YfKdwxpOrieUIIIdxv//EivtmYpinzMum5Z3i7OvkRmXl8B8u2fEhsbiqD8K62TjvFrHluV1WoIdkwGMzYnfZqtzWkZpuIRIV0pF9+V46VpJNRnks0Nd9+2aSWYEcl317CorTlLE1bwctjPsI7on6SEUVRmDogjiN5pRw5pXlvQ3IubcP8qoywEUII4V75pRW8v/JglZGOtw5qTUTAhd1WTz+2ifkb/49t+ftRUQnGwAC9F/pqfpS2wIg3OqxAS+8wIsxBoFTfHTQsuAN9W4+5oNjqQrNNRBI6XkZCx8sAcNgrsBakQXEm5Kf+9TgMJZUznK5XtZOLdTGF4B3WqV7j8zJWZtH//GU3Zad0ePp2YxpxIT50aIAOT0IIIc7O7nDywapDFJy2hsyErlH0aBl0wccvLs1ia/4+1/Nc7GxRS+mt+J6spDNAaHt04Z24zy+E8BZ9MZnP/D0RFNyWXsFtLzi+C9VsE5FT6Q0mfELaQkhbiBtwckN5AeWHVmDb+TGUn7zn1z9uOOjqf8BReIAXtw5uwzvL9rvKnE6VD1Ye5JmJnWQImBBCeIDvNx1h//EiTVmn6AAmdW9RJ8dv12Y0cVve53BZpqtsJSX0juwP4QkQlgAh7cBQ+Z0QUydnbTjNevjuWXlZ8Oo0iSeumc9jw99kWNQAogy+xCdcXeMuS39/iSW//wubrbROQugeG8gl3aI0ZQVlNj75LRm1hiFZQgghGsa2tHyW7D6uKQv2NXH7kDbodLXsz+d0Vlus6HRcHH8NAF46I6NaDOG2S2bCiH9A5yshopMrCWmMpEXkHMXGDiA2dgCqw46ir/7Plpd7kPnJ87GpDtakrWRSl5vo1ulalAtsPbmsWwuSs0vZdbTAVbYnvZDV+7MZ2iHsgo4thBDi/JRW2Pl8bYqmzKBXuHt4O/y9zn3ah5LiDL5d+RQ9AtrQfchT1dbp2ularinNpneX6/HxCb2QsD2OtIjUUk1JCMBPf7yETa3sz5FtK+LTzf9H9tKn/hpadf50OoXbh7QhyFeb8X67IY3ckooa9hJCCFGfvtmQVmVF3Sl9W9I61LeGParases7XvzpOjbl7mZOygKK0v6stp5Ob2BIv/uaXBICkojUmazMXWzK3a0pG6T4EZZ9CBY8Bn9+DOUFNex9dn5mA9MGtNKUldscfP5HityiEUKIBrbzaAG/78/WlCW2sJx7K3VFCYdWvsiHm96k0FE5OrIYB9/88TxqRd3c2m8sJBGpI2HhiTw89BXa+EQC4I2O8TrLX1tVOLAU5j8E6dvP+xxdYiwMbKfNhnceLWDtwZwa9hBCCFHXTvwIPJXZqGPagLhzm+cpfRv8+gitj26ji6KdD2SfrZC8rN017Ng0SSJSh+LiBvPAFXOZ3vUOrvRpib9y2gQ21iJY8S/Y8X2Nc/+fzeQ+sVXWo5n9Zyr5pXKLRgghGsJ3m45UuS1+de9YQvzMNezxF1tZZev4in9BaQ6KojBFF4zPX1/FiQFteOqSzwlu0bu+QvdIkojUMUWno3f36fSfNAu6XA3603syqxRsn8PnP0+juCi9ukOcka/ZwNT+cZqysgoHX6w9LLdohBCinu3NKGTl3kxNWXyUP8POdksmKwl+fbiydfwUAYqe64yRXJ84jTsv/QJLYKs6jtjzSSJSXwxm6HIVXPJvaNHLVexQVWY6s9lQsJ9X/jeVlJRVtT50j5ZB9G0drCnbmpbPn8m5NewhhBDiQpXbHMxak6IpMxl0Z11HZt+u77AvfQ5KsqtujOxC90s/ZECvOy94hGVj1TyvuiH5hsCQR6D7dYDCz2oBB1QrAHn2Et767Ul2bvu81oe9rl9L/L20I3i+Wp9KYbn7FzASQoimaN6Wo2QVWTVlV/aMIdy/5incV6x9nXc2vcU39kxtq7XeBL1vgeFPgW/TGwlTG5KINARFgU6XUTr0UTYo2vuKAehoHdal1of09zJyw2m3aEqsdr5al1rDHkIIIc7Xgcwilu3RTlzWLtyPkQnh1dZXnU7mLX2EufvnoqKyVi3hF/WvkZNh8TD+NegwuvL7oZmTRKQB+bToyWMTZtLONxoAPQq3JN6Mb3TP8zpe71bB9IzTrmOwMSWXTYflFo0QQtSVCruTz9akaMYYGPU6brqodY23ZDKzdrI6fa2mbJGzkP2xPWDkDPCPrM+QGxVJRBqYJbAVf7vsa0a1GMKVEf1o1WP6BR3vhv5x+Jq1t2i+XJdKsdX9SzsLIURT8NPWoxwvKNeUTerRgkhLzbdkIiK6clPPv6FwMlGZ2Ho87S56tEHWKmtM5K/hBnqDiUkjX2HIxW9ccLOcxdvIlL4tNWWFZTbm/Cm3aIQQ4kIdyipm0a4MTVnrUF9Gd4o4675dE6/l2vgp6NBxQ+J0xgx+utl2SD0TWWvGnc7whkza9wtWWwldE68962H6twlmQ0ou29LyXWVrD+bQp1Uw3WID6yBQIYRofmwOJ5+tSdbcktHrFG4e1PqcF7Qb1PdvdGw1grDwxHqKsvGT1MwDHTi4mA/Wv8wnm/6Pzdv+e9b6iqJw44A4vE3aCdT+u/Yw5TZHfYUphBBN2sKdGaTna2/JXNo9muhA7WyotuqG5Z5CkpAzk0TEwySnrOS9P57Hpjpw4mTmtg9Yv/mjs+4X6GNich/tLZr80goW7syoYQ8hhBA1yS2p4Ncd2kknY4N9GJuo7WSanLKS5+ddyfGkXxoyvCZFEhEPs3n//6hQT3Y0VVHZvetb1H1LzrrvRe1C6BQdoClbuDOD7GJrDXsIIYSoztxNR6iwO13PFQVuvqg1Bv3Jr82MjK28//sMcp0VvPfnqxSkra3uUOIsJBHxMFeMfI0hkf1cz7sq3kzVBaNs/ASSFpxxX0VRmNy3pWY4mc3h5LuNR+otXiGEaGoOZBax7pB2MdEhHcJoGeLjel6Qn8J7Kx6m1Fk5N1SOauP9VU9RnpfcoLE2BZKIeBhFp+PqUW8yssVgOile3KQLxXAisdg0C3b/dMb9WwR6Mzxeu+bBxpRckjKK6iliIYRoOlRV5es/0zRl3iY9k3q0OFlQUYL1tzfApm1tDvCNwOAX1RBhNimSiHggRadj0vCXub3rHRhPH967dfZZV++9rHuLKnOLfP1nKk6nLIonhBBn8sfBHFKySzRll3aLJsDrr1XP7RWw+jXCizJ5SB9BtFJZ3tI7nFvGfYDBWPPcIqJ6koh4KEWnw9BtMnSbXHXjju9g17wa9/UzG7ise7SmLC23lN8PnLlntxBCNGflNgdzN2lvZUdYvBgR/9c07k4nrH0HMvcAYFH03K+LoL93FHeNeQ+zl6WhQ24SJBHxdImXQ88bqxSv3PoJG7Z+VuNuwzqGVxli9sPmI5RWyIyrQghRnV+2p1NQpl04dEqflpUdVFUVNn0GaX9qtvt4B3PDuA/wD2iBOD+SiDQG8RMqV2kEnKrKd848vnfm8eWOTzl4sPrRNHqdwuS+sZqyonI7v2xLr7a+EEI0Z5lF5Sw+bQbVLjEWusT81cqx6wfYf9rnrcELhj0OftUvfCfOjSQijUWH0dh73MCHzmxWOSs7njpUJx+ufYHM4zur3SUx2kL302ZWXbrnOBmnrZkghBDN3Xcbj+A4pR+dTqe45mZat/F9dmz9XLuDzgBDHoHg1g0ZZpMkiUgjYkiYSGhEN01ZqbOCvTu+rHGfa/vEoj9lKmKHU+WbDWk11hdCiOZm97FCNh/O05SNjA8n0uLFjl3f8dXuL/jImcUfzuK/tiow4F6I7NzwwTZB9ZqIvPTSS/Tp0wd/f3/Cw8OZNGkSSUlJ9XnKJu/Kka+RGNAGAD0K01qMZMjwF2qsHx7gxcWnLc60/Ug+O48W1GucQgjRGDicKnM2aBcJ9fMycGn3aA4f/o3PNr+NiooKzHbmstBZAL2mQdwA9wTcBNVrIrJq1Sruuece1q1bx5IlS7DZbIwePZqSkpKz7yyqpdMbuHnse3T0juTvCTfSZ8TzoD/z2oWXdI0mwNuoKZuzIRW7w1nDHkII0Tys3pfF0bwyTdnlPVrgYzKwdPtn2FTtel0V0T2h47iGDLHJU1T1DBNS1LGsrCzCw8NZtWoVQ4YMOWv9wsJCLBYLBQUFBAQEnLV+s+Kwgd549np/+W1/FrPWpGjKpvRtyahzWMpaCCGaomKrnSd/2EGJ9eRowpggb2ZMTESnU7Dbyvlp5VOsSP8DgH6hXblh7PsoZ1g5XVSqzfd3g/41CwoqbwcEBwdXu91qtVJYWKh5iBrUIgkBGNQulLgQX03Zj1uPUlRuq2EPIYRo2n7eekyThABM6dcS3V/96gxGL668+A1u7fk3EgPaMOXityUJqQcN9hd1Op3cf//9XHTRRXTuXH0Hn5deegmLxeJ6xMbGVltPnFlpaTapqb9ryhRFYcppw3nLKhz8uPVYQ4YmhBAe4Vh+Gcv3ZmrKesYFER9Z9dd7987Xcdekr2TW1HrSYInIPffcw86dO5kzZ06NdZ544gkKCgpcj7Q0Gd1RWznZ+3jzfzfyzuonycjYqtnWPsKfvq21rVGrkjI5lq+9PyqEEE3ddxuPcGrPBINe4Zre8uPXHRokEbn33nv55ZdfWLFiBTExMTXWM5vNBAQEaB7i3B0+/BuvL7qDDGseZU4bH618nPIy7ZC0q3rFYDxlGWtVrbxFI4QQzcWBzCK2H8nXlI1OCCfMW6l+B1Gv6jURUVWVe++9l3nz5rF8+XJat5aJX+rT8h2zKHKcnKwss6KA71c8rqkT4mdmdKK2g+qmlDwO58hIJiFE06eqKnM3a398+XsZ8Cn9gq9+vA5bXrKbImu+6jURueeee/jyyy+ZPXs2/v7+ZGRkkJGRQVmZ3AqoD9eNepNYr1DX81AMDM47DodWaeqN7RyJt0mvKfths7SKCCGavt3phezLKNKUDYo8wK8p81lbepTX5t/M8b0/uym65qleE5H333+fgoIChg0bRlRUlOvxzTff1Odpmy2zl4Wbh72MWWegl+LDY/pI4hQzbPwUCk92SvUxGRjXOUqz786jBRzILDr9kEII0WSoqlrlR1eouYDNRz/4a8oyOOYs5/U/X6Eke587QmyW6v3WTHWP6dOn1+dpm7Ww8ESe6PsE03UheCt/vbx2K6x5u3Lukb+MTAjH30s7EdrczUdpwGllhBCiQW1Nyycl+5Tb0KrKWN2PlDu0rfTD4kbhG9qhgaNrvmRAdBMU2mE8SruR2sK8FNhyck0aL6OeS7pGa6rsyyhi1zGZu0UI0fQ4nSrztmhbQ0aqfzDUnsmj+kiilMq5mTr4t2TcoKfdEWKzJYlIU9VrOgS00JbtWwhHNrqeDu0YRpCvSVNl3hZpFRFCND1/puRqpnJvUZHMWOdqdIpCpGLkEV0Ew00RTB/xOrqzLJsh6pYkIk2VwQwX3VdlBtb9f7xBXu5BAIx6HZd207aKpGSXsCUtv6GiFEKIemd3OPnplGkKDM4KLiv/iSCfkwmHSdFz5fB/EWCRuUQamiQiTVlQHPScDoBTVVngLOD/ylOYteJRHPYKAAa2DSE8QDtb4I9bjuJ0SquIEKJp+ONgDpmFVtfzASVLaedbhsIp84Z0vQYiOrkhOiGJSFPXbiSFLXrwjjOT+c4CVOBgyTEW/P48AAa9jkndta0iR/PKWJ+c64ZghRCiblXYnfy87eSowQhbGgPtfxLgfcrtl7B4SLzcDdEJkESk6VMUDH1uJcegvee5KHUZhw8tBaBv62Bigrw123/edhS7w9lgYQohRH1YtS+LvJLKFmBFLcdpe59SCydbQ/RG6HcHKDKrqrtIItIM+PiEctOAf6D76+VWgEuCuxAb2aPyuaIwqYe2Y2tmoZU1B3MaOlQhhKgz5TYH87efbA2xWL/kuDGfjw25zHHkUqY6ocvVEBB9hqOI+iaJSDPRutUwJraZQJBi4IFONzFmwgfofEJc27vHBtI61Fezz/+2HaPCLq0iQojGaeme4xSV2wHwtu0hW9mM2ahHAX5Xi/nc7IT4S9wbpJBEpDkZNfBxnhj7CW16316lGVJRFC7vqW0VySupYGWSdplsIYRoDEor7CzcmeF6bnIsRK8Hg67ys0+Hjkv6PwI6fU2HEA1EEpFmRNHp8AnrWOP2TlEBxEf5a8p+3ZFOuc1R36EJIUSdWrgzg7KKk59dBea/cVnbKzEplYnHxbHDiInp767wxCkkEREuiqJweY8YTVlRuZ2le467KSIhhKi9wnJblc+t3q3DuWz44zw57hMuCu/F2IFPuik6cTpJRAQANlspGelbaBfuR9eYQM22hTszKLHa3ROYEELU0q/b07HaTvZvUxSY1KOyQ2poaDxTxv4Ho9m3pt1FA5NERHAoeRmvzL2Sd1Y8TGlpDpefNoKmrMIhrSJCiEahoNTGitP6tg1oG0qUxbuGPYS7SSLSjNlt5Xy76D7e+u1pMiryKbCX8r/fnqNliA+9WwVr6i7dkyl9RYQQHm/R7gzsjpMzQ+t1SpWlLIRnkUSkGdPrTWSWHkfl5D/a345v4MDBJUzsFqWpW2q1s2KvjKARQniuYqudlUmZqKqTFmUfE1axn4FtQwjzN7s7NHEGkog0Y4pOx5TBz7l6kQPEKia8d3xPjL+BHi0DNfUX7cqQeUWEEB5r6e7jWG1OAm0LSdZtJt/5b9qX/Rds5e4OTZyBJCLNXEhoRy5peylGFC7XBfKwLoIWpfmwax4TumqbM4vK7azel+WeQIUQ4gxO9GXTOwspVBcBUGGAr48tZM4vN4MqC3l6KklEBMP6Pcg/QgcwUheA/sREZ7t/orU+m8ToAE3dhbsyZA0aIYTHWZGUSVmFA/+Kr7FSubaM2Vj5Fdex7XhZS8aDSSIi0OkNhAz4O5xyiwbVAX9+xCVdIzV180oqWHtI1qARQngOq93B4l0ZBDjyGF2cSku7EaNeh15RSAxoQ/fO17k7RHEGkoiISkFx0OlSbVnOATrk/077iKqzrTqc0swphPAMv+3LpqjczkXFi4hwKFxV7M+dhlAiFRNXDHgMRSdfdZ5MXh1xUuIV4K8dLcO2OVzSVtukmVloZUNKbgMGJoQQ1bM5nCzclUF0RQrtyncCEOBlpK/Jn6cSphMR0dXNEYqzkUREnGQwQb87XE9LVSdfV2SwctcMWgZpJwOavz0dVTp/CSHc7I+DOeQVWxlc/KurLCLAC0x+KF2ucmNk4lxJIiK0whOg3Si2OUt50ZHOGrWYvUWH6WpZral2LL+MLWn57olRCCEAh1NlwY50Esq3EG47BoCf2YCf2QBdrgaz/1mOIDyBJCKiioouV/GtrpQCTs6kuunYD7Q47d/0L9ukVUQI4T5/JueSW1DIgOIlrrLwADMEREO7UW6MTNSGJCKiCpN3EFd3uVlTVmgvpb/lT03Z4ZwSdh0rbMjQhBACAFVVmb/jGP4V/2W1VzpFigMfkx5/LwP0vBH0BneHKM6RJCKiWt0Sp9Det3Lxu0RTEE/1eYwxw++q/LVxil+2p7sjPCFEM7c5NY+cnAPkKFvZbbLymaWA7RYbtsguEN3D3eGJWpCUUVRL0em4ut8jFBxZR0KvO8DoBcD4LlHMWpPiqrf/eBFJGUV0jJR7sUKIhqGqKr9sT8ds+5JipXKCRVUHqwwlDEqYQORZ9heeRVpERI2iY/qR0P8+VxICMKBNCEG+Jk29+duPNXRoQohmbOfRQtKzDlOoHHWVmQ06Bkf0JjKqpxsjE+dDEhFRKwa9jrGJ2t8bu44Vkpxd4qaIhBDNSWVryDFsujBMxmcIUzugVxQCjGbGD3zC3eGJ8yCJiKi1IR3CKjuEnUJaRYQQDWHf8WIOZBYDYNOFUeD1ANd0e5Hru9+J3+kTMopGQRIRUWsmg46BEQdQ1DJX2ZbUfI7klboxKiFEc/DLaT96gnxNjO49StaTacQkERG1kp29l0/+N51VqW8S6vhes23Bjgw3RSWEaA5SskvYfdqUAeO7RGLQy1dZYyavnjhnf2x8j+cX3MLWvCQUBQp0GzE6jru2/5mSS25JhRsjFEI0ZYt2aX/sBHgbGdQuzE3RiLoiiYg4ZzHhXXGeMpOqTu/E1/6V67nTqbJ09/HqdhVCiAuSXWwlJWklRqfVVTYqIQKTQb7GGjt5BcU5a9lyEP3CTq5k6avouUjNItR28p7tqn1ZlFbY3RGeEKIJW7BpM+nqhwSWP0mb0p/x1dsZ1lFaQ5oCSURErUwc8Dg+OgPDdf48q49mgk8gQ4t/hb9aSsptDlbvy3JzlEKIpqTEamd70gc4cHJUX8pe/QKiKmbgo3OcfWfh8SQREbViCWzFc13v4UpdED6KDrNBTyfDUdpad7vqLNmdid3hdGOUQoimZNHmjWSpe1zPFSAqoh2K0VzzTqLRkERE1Jp3pyvAN9T1PMzfzKDihejVylsy+aUV/Jmc667whBBNiM3hJH37N3ifstC3WW/g0oEPui8oUackERG1ZzBB9xtcT31NBqINRXQrXesqW7QrA/WUjq1CCHE+Nu3ex4CCJG4pCGRYmQ++qo4B4T0JDm7n7tBEHZFERJyflv0htIPraZi/mV6lq1092o/klbHrtPH+QghRG6qqkrVuDjrVgQGFXlZv/lYRyeShMpV7UyKJiDg/igI9p7me+nvpyTAXElJxcpKzhTtlgjMhxPnbvf8AMbnrNWWWTmPxscS6KSJRHyQREecvtB1qdC92OMt41ZnJr34lpCvr0DvzAdiTXkhqjkz7LoQ4P+lrZqNwsuO70WSmxcDJboxI1Id6TURWr17NxIkTiY6ORlEUfvzxx/o8nXCD4x2G86EziyNqBUaDDofiwP/UVpFd6W6MTgjRWKUePkRI5lpNmTl+DIpPsJsiEvWlXhORkpISunXrxrvvvlufpxFuFBndh+5BHYHKIXUmg45cZSsGZw4AfybnkVNsPcMRhBCiqnWr38V5SmuIzmAibogsbNcUGc5e5fyNGzeOcePG1ecphAe4pPd9bFtyDyoqMUYvzNaupOKFncrOZkt2H2dy35buDlMI0UgkJW/ml7LVBATo6F/uTacKE7oOYzD6Brk7NFEP6jUREc1DZFQPLo7oQ6TJQu8+f2POrjL27zm55szq/Vlc2j0aH5O83YQQZzd3zZuoqkqBzsEin2I2e5l4fsgUd4cl6olHfTNYrVas1pPN+IWFMvyzsbh09L8rR9IAF3eysnzv8ROzvmO1OVmZlMX4LlHuC1AI0SgcSt1IUmmSpiwmuBe+lhA3RSTqm0eNmnnppZewWCyuR2ysDNFqNP5KQuCvOUXitB3Klu45LtO+CyHOamNyBv7Ok7dgDOi5atgjboxI1DePSkSeeOIJCgoKXI+0tDR3hyTO09jOkZrnBaU21su070KIM6iwO1mbEUOp+XlClOsIUIOI9+9DywjpY9aUedStGbPZjNksixg1Ba1DfekQ6c++jCJX2cKdGQxsG4JySuuJEEKcsD45h8IyG4qio8g0GFW9iKtHylTuTV29togUFxezdetWtm7dCkBycjJbt24lNTW1Pk8rPMTIdt4EWmcRWP4eAMfyy9h5VPr9CCGqUlW1ymzMCdEW2kRI35Cmrl5bRDZu3Mjw4cNdzx98sHK1xGnTpjFr1qz6PLVwowprEas3/IfFKQvJ0ZXidIK/PZlyQ2sW7kqnS4zF3SEKITzM9iMFZBSUa8rGJkoH9+agXhORYcOGyQqszVBZeR7zk+djUx2YDDrKKhyY7N9TbniEvelFpOWWEhvs4+4whRAeZMnu45rnLYK86dwiwE3RiIbkUZ1VRdNgsbRkSNQAAIx6HTpFIVdJxtd+GKgcQSOEECesW/c52Yc3wSk/XC/uFCH9yZoJSUREvRjV935MigEF6Kp4c31RAINK1gGw7lAOheU29wYohPAIpaXZzNr1IVnOt2hV+gIx5RsIMOvo11r6hjQXkoiIeuEf0ILJrSfwsD6Cv3tFEOE00KF8ByG2DOwOlZVJWe4OUQjhAZb/+QGljgoqFJV9hmOkKjMZEluBySBfT82FvNKi3vTtcy+tTEEY9TqCfEwA9CtZDsDKvZnYZIIzIZo1u62cpclLTr0jQ6Aaw4ievd0XlGhwkoiI+mP2g/gJAIT5Vc4P09a6h0B7NgVlNjakyARnQjRnxQeXE15q59SeIB1jr8XiY3RbTKLhSSIi6lfH8WD0wdukx89sAFR6lv4GwOJdx2VUlRDNlapStuVXxhf7cXNhIN2tXoQ6wxjf/1J3RyYamCQion6ZfKD9aKByDRqAhPKt+DoKScstZX9msTujE0K4iXpsK4XphwAIdOoZWeZL99BbaBkiQ/ubG0lERP3rOA50BgK8DJQaVZZ7FWKu+BqoOneAEKJ5yN00lzKbw/U8zxBKl74j3BiRcBdJRET98w4kK7YPnzlzmBVYwDZzOdm6Heid+WxJzSOryOruCIUQDSnnIIUpWzVFB0NG0D02qPr6okmTREQ0jPYXs1Utw6jXoSgKdhz4V/yMqsIymeBMiGaleOs8CstOziVUpvMjrtdodDqZwKw5kkRENIiw8ES6B3VEUcCkr/ywKWYzOtXKb/uzKatwnOUIQoimoCD3IAf2reDUbuq7/QcyqGO022IS7iWJiGgwF3e7BQA/g4H29mjClRtxYqLc5mDNgWw3RyeEaAgLNrzDm/pMfvQt4ojBhk0xYuk6AW+T3t2hCTep10XvhDhVy5aDuLHNJDq3Gc3MvX5sTs13bVu29zgj4sOlaVaIJqykOIPVx/5EVVUOGis4aKwgztmWv3dp5e7QhBtJi4hoUH0HPYZPdA9GdYrUlGcWWtl2JN89QQkhGsTqzR9RYj/ZN0RBISLuGsL9vdwYlXA3SUSEW3SI8CM2WDtfgKzKK0TTllviRO88eQsmWG3J6J793RiR8ASSiAi3UBSF0Z0iNGV704tIyy11U0RCiPp2TDcZvfFfhDMYb9WLkKAr6BDh5+6whJtJIiLcpk/rYCze2jUlFssEZ0I0SWm5pexNL8Kp8yfffB0O82uM7DUGRZF+Yc2dJCLCbYx6HUM7hOJf8Rs+5c/gY9vB+kM5FJwyv4AQomk4/UdGgI83fVoHuyka4UkkERFus2PXd+xMupNcZpOvZKF3/ILDqbIyKdPdoQkh6lBBmY31h3I0ZcPjwzHq5StISCIi3Ohozh7SK3JcH0a5Shpe9mRW7M3E5nC6OTohRF1ZmZSJw3lyCjO9TmFYxzA3RiQ8iSQiwm2G9LoLk2LA9FcioqLia/sfReV2NqTkujk6IURdyEz5nbW7UjRl/duEEOBlrH4H0exIIiLcxsc3jMFR/dHrFMIUAxeX+nJlUQ7ejmKW7clEVdWzH0QI4bGs5QW8sPJxKgofpk3ph4TYkgEYlRBxlj1FcyKJiHCrEb3v4RZDBE/ro+ha4YVZddC9bC0p2SUcyi5xd3hCiAuwcfsXFNmsZOls7NVvJd/xOp2Cy2gZ4nP2nUWzIYmIcCtLYCt6tBlLoJcRk6Hy7di1dD0GZwVLZSivEI2W6nSy9MB87Kf0DfEmgmHdOrsxKuGJJBER7pcwEUXREepnBsCklhNfvpWNh/PIL61wc3BCiPNRlLGVwpJ8TZmf90i6xwa5JyDhsSQREe5naQHRPQj2NXFizbtuZetwOpysTMpyb2xCiPNi3LeKG/L8ubzYn9Y2E76qDwO7X41eFrYUp5FERHiGDmMw6BSCfE0ABNszibElszJJhvIK0eiUF5Cf9BuKqtDGbuKKEn+6qJcxNF46qYqqJBERniGqG/hHEepnxo7KDlM5DvsPMpRXiEbIsW8puUUn141yKAYCOo3Bz2xwY1TCU0kiIjyDolDSdiiL9EV8FlzAYp8SUvWHMDtSWbpbhvIK0Wg4neRsm4/NcfLf7D5zF4Z1aePGoIQnk0REeAx9q6H8RikOQ+U9ZBUVb9uvHM4p4WCWDOUVolE4tpn87AxNUVHLEcQGy5BdUT1JRITH8PIJpn94bwx6Bd1fHdryld0oahnL9shQXiEag+wtP1Na4XA9zzS2oFfP3m6MSHg6SUSERxnS7SYUwKTXEePwobVjIAZVlaG8QjQCmcd38nzGYtZ4lVKkVCYjyUGDZMiuOCPpOSQ8SkREV65pOZq4kK68trUFmQ595QanysqkLCb1aOHeAIUQNVq+5TPyVDvrvOz86VVOB5svvXqOlCG74oykRUR4nCHDniOuy5X0bRepKZehvEJ4rgprEb8dW8+JfuVOVPKNbRgULz8exJlJIiI81sjTFsYqKrezIVmG8grhiXbtm0+xXXv7tGObyTJkV5yVJCLCY7UI9CY+yl9TtlRW5RXCI5X7XIy/8jfC1HgM6AlSIxnTa5C7wxKNgKSqwqONSohgb3qR63nlUN5i2oX7n2EvIURDW7bnOGXGBMqMCRidhcQEV8iQXXFOpEVEeLRuMYGE+JkwOY4SWP4ZilrG0j2Z7g5LCHGKg1nFJGefnOvHoQtgZI/+boxINCbSIiI8WnLyUoIrZpJh30+xohJeEcKmw5PIK6lwrUsjhHCv0+f5CfEz0T020D3BiEZHWkSER5u/7WPSbMmgVPYLKVfX4nA4WLlPWkWE8AT5pRVsTMnTlA3vGC5DdsU5k0REeLShHa9EUcCor3yrFioF+Nm3sCopiwq7DOUVwq1UlVVJmTicJzuQG/U6BncIc2NQorGRRER4tC7xVxJs9MVk0GFEoWuFF93LdsmqvEJ4gJ175rFt3S20KpuPyVm52u6AtiEyZFfUirxbhEfT6Q2Mbz2BkgOLiSgHe7kTlQNs9S1k2Z5MBrYNQVGkCVgId5i/fQ77dTnAL4SVLyZCGcSI+GfdHZZoZBqkReTdd9+lVatWeHl50a9fP/7888+GOK1oIvr3uJWRxhBi/bwBUHDSuWyDayivEKLhZWfvZW9Riut5lr4Cs8VHhuyKWqv3ROSbb77hwQcfZMaMGWzevJlu3boxZswYMjOls6E4R2Z/aDWYAC8DJkPlW7ZL2QZ0qoNlMpRXCLdY8OfH2J0n+2kZMTC8z01ujEg0VvWeiLz55pvcdttt3HTTTXTq1IkPPvgAHx8fPvvss/o+tWhKOoxFURRC/cwAeDuLaV++Q1blFcId7BV0PZzK4DIfApyVC1NG6OPp07almwMTjVG9JiIVFRVs2rSJUaNGnTyhTseoUaNYu3ZtlfpWq5XCwkLNQwgAguIgPIFgXxMnRgV2K1uP869VeYUQDad4/+9QVEhfqze3FFq4vNiffgnTZMiuOC/1mohkZ2fjcDiIiNAuXhYREUFGRkaV+i+99BIWi8X1iI2Nrc/wRGPTYRwGneKayMzpTMbftltW5RWigWVs+p9rlV0dCt76BEb1HuzeoESj5VHDd5944gkKCgpcj7S0NHeHJDxJTG9s3kEk+9r5yr+A2f4FKI5fZCivEA3Inn0Ia8YeTZmz/cUyZFect3p954SGhqLX6zl+XDv97/Hjx4mMjKxS32w2Yzab6zMk0Zjp9KwOi2Ve0TZKTA5wQK5yGG9nFsv2+DKgjQzlFaK+HfnzJ2yOkxOYlej86dpvpBsjEo1dvbaImEwmevXqxbJly1xlTqeTZcuWMWDAgPo8tWii+nW7GYOiw6Sv7CDnxImv7RdSsks4dMqiW0KIelBRQtn+VZqinKghxIYGuCkg0RTU+62ZBx98kI8//pjPP/+cPXv2cNddd1FSUsJNN8kwL1F7fv5R9ArpgkGvoFMUApx6LA47UHXhLSFE3Vq+9hO2O/JxUNkioqIjts8lbo5KNHb1flPv2muvJSsri2eeeYaMjAy6d+/OwoULq3RgFeJcDesyjbzfZtAitC+r8wZxWB8CwMaUPK7pXUGgj6zKK0RdU51O/nfwJ7J9i/FVdXSxmgk39uKK9q3dHZpo5BRVVdWzV3OPwsJCLBYLBQUFBARI0584hb2CYoeOh7/dphkxc2n3aC7r3sKNgQnRNG3d/TNv/v5PTv3GGN/2AaaMut59QQmPVZvvb48aNSPEOTOY8DMb6N8mWFO8MikLuwzlFaLO/bpltiYJ8cOPCQOvcl9AosmQREQ0aiMTtLf4CstsbEjJc1M0QjRNdoeTEltXAtUQV1nHoCEE+MgoR3HhJBERjVpssA8dIv01Zcv3SqdVIerS5tR8MriYUq8XsOj+TpiawPiBt7o7LNFESCIiGr1RCeEAKGoF3vb9HMoq4ZCsyitEnTl1RFqZMYGwlk/SIUbWlRF1QxIR0ei1CSgmwjELxfooFY73QLWxfK+syitEXTicU8KBTG1if/otUSEuhCQiolHLzt7Ls79MIUvZSLlipYxyAmwr+DM5l4JSm7vDE6LRW7pHm9RbfIz0bBnonmBEkySJiGjUQkPjifEKxaRXODG7u935Gw6nyqr9siqvEBeisNzGn8k5mrLhHcMx6OWrQ9QdeTeJRm9ou4koioLxrw9HP0cBgbYjrEzKlKG8QlyARcs/ILR0JyfG7ep1CkM7hrk5KtHUyHKJotHr0fk6ftk9mzZOheg8CHcY2OW9keXGGDYdzqNfm5CzH0QIoVFUlMnCI1/hxEHL0iAMuotoEX8lAV5Gd4cmmhhpERGNntHowzOdb2e6KZTWhsp5DeLLt+LlLGWprD8jxHn58bf3qFDt2FE5ZMhlv+4XesdKC6Ooe5KIiCbB0H40KHrC/CoTEb1qJ6FsswzlFeI82G3lrDu2XFPWwhhH13Zd3RSRaMokERFNg08wxPbF4m3E9FdfkW5l61BUp7SKCFFLqTuW0arMiQHFVTa0041ujEg0ZZKIiKaj4zgURSHEr3L1XX9HPq2te9mQkkd+aYWbgxOikVBVSjf/j4tLfbmtIJAB5d60IIJRfS5xd2SiiZJERDQdoR0gqBUhfmZ0f/2Qiy9bidOpsjJJhvIKcS6KDm/DkZsMgI+qY2C5DxM73C1DdkW9kXeWaDoUBeInYNApVPgqzPcp5mefnfjaNrIyKZMKu3S0E+Jsjq77VrPKbpEhhK79R7kvINHkSSIimpaWA/lKb+VTv3z2mqw4UVEcv1JUbmdDSq67oxPCo9lyD2NP26wpK287Fn9vWWVX1B9JRETTojcQ3XIQekXBoK+8P5OrpONj38uS3cdRT/2pJ4TQSPvjO+zOk/9GynXexA+c6MaIRHMgiYhocgb2uA0fnQmTXg+AEYVg22bSckvZnylDeYWoTmFBGsmHlmrKsiOHEhse5KaIRHMhiYhocsxeFobEDCFQr6cL7YhU/sYRrykALNktQ3mFqM7Pa/6PWeZsvvYr4ICxAruiI7b/5e4OSzQDMsW7aJJG9bmP0YlTWZkZwDcb0lzlW1LzyC62Euon97yFOMFaXsCa9D8AOGaw85OhiDa6lsxoG+fmyERzIC0iokny8g3FFNaBQe1DMRtPvs1VFVbszTzDnkI0Pys2fEyJ3aop695pKjqdUsMeQtQdSUREk+ZjMnBRu1BN2er92ZTbHG6KSAjPk5rrxKyebCUMIYoxfYa5LyDRrEgiIpq8kfERmuelVjvrDuW4KRohPIvV7mBn6cWo5lcJYwL+qoVOsVfhY5I796JhSCIimrxIixedW1gAUFUnRsdxlu6RobxCAKw7lEuJ1Y6qmCgwX0K5+V9cMuAqd4clmhFJRESzMDI+hICKxXhbn8Jpe4X0vCJ2pxe6Oywh3EpVVZadtihk55ggooN83RSRaI6k7U00eUWFR5n3+83kKbk4lMpWkDDbYpbuDiMx2uLm6IRwn70ZRRzNK9OUXZwQUUNtIeqHtIiIJs/PLwofgxemUxbtKneuZvuRXDILy90YmRDutfS0eXUiLF50bhHgpmhEcyWJiGjyFJ2Oi+OvwajXoSiVwxH1agmR1h0sk6G8opnavuVbklOSNGUj48Nd/0aEaCiSiIhmoWuna4kwW4jTmbikxI+bCi0MLN3I7/uzKauQobyieSkuOsa7m9+ivOJ5Wpe8RVjFfryMuipD3YVoCJKIiGZBpzdwX/d7+IcxinibGR0KLSqSsZQeZs2BbHeHJ0SDWrT+bcocdop0DpIM+8hx/pvBLcHLqHd3aKIZkkRENBuW9mMxB4QS4G10lfUqWc3SPcdxOmUor2geSorSWZ72m2b4uj+tGNE90Y1RieZMEhHRfOiN0HE8Yf4nZ5Bsa92DLe8oW9Ly3ReXEA2oZM//CDulj7aCQqsW0wn393JfUKJZk0RENC/tRuHrF4C36UQTtErvklUs3p3h1rCEaBDWIoy7lzOp0J/riiy0thsJUVsxps9gd0cmmjFJRETzYvJBaX8xYaesvtvaupn0o7s5lFXsxsCEaAB755OdXzmRX5TDwBXFFmJC76RDhJ+bAxPNmSQiovmJn0Cgvx+qXmGdVxmfBOSht3/NktPmVBCiSbEWUbz9f5RYT44S2+fVhYt6dpchu8KtJBERzY+XhYNxvZkZlM8ar1KsikqOcpidB9aRU2w9+/5CNEZJC8jOLzilQCEp9GJ6xwW5LSQhQBIR0UxFd7kWg8nIiR+CKip621yW7ZEJzkQTZC3GuusXCkptrqJ9Xp3p3bULBr18DQj3knegaJZ8/SIZETsM4ykfwv6OAtbsTaXcJhOciaZl5+ZPOJiXw8kBuwrbLCMZ0iHMjVEJUUkSEdFsjeh7HxaTmZYOfzo4ryDD50WKHSZ+2y8TnImmo7Q0m5kHfuA1QyZLfEoo0DnY79WZxIREfM2y7qlwP0lERLPl4xPKc6P+Q4sO75LifTFOpXJI75LdGTLBmWgyVm34D4U2K3ZVZbupnJkBBfzh149RCeHuDk0IQBIR0cz5R3VjVJdYTVlOcQVb0vLcFJEQdae0NJtlacupcDhdZYFqLO3a9CA8QCYwE55BEhHR7LUP96NVqK+mbPEuGcorGj/VaaedXxcU58nhuXb9JEZ3inBjVEJo1Vsi8uKLLzJw4EB8fHwIDAysr9MIccEURanywXwgs5iDMsGZaOR8/SLJMd6Dt+EfhKmdCVVbExnVm3bhMoGZ8Bz1lohUVFRw9dVXc9ddd9XXKYSoM73iggjyNbmee9n3sWD9YjdGJMSFO5BZzKGsEir0URR43UOB+WFGJ0bKBGbCo9Rbl+nnnnsOgFmzZtXXKYSoMwa9jlEJEfxv3UpM9rnkKCkkHQkmq3AcYQHe7g5PiPNy+hpKIX5e9JIJzISHkT4iQvwlWr+BYucbZCvJqKjkKzn87/ev3R2WEOclq8jK5sPaTtcjEyLQ66Q1RHgWj0pErFYrhYWFmocQDSWh3UiCjNpOq1vSvqa03FbDHkJ4rmV7jqOeMgrdbNQxpEOo+wISoga1SkQef/xxFEU542Pv3r3nHcxLL72ExWJxPWJjY8++kxB1xGj0YVyHKzjxe9GIQnx5OVv+XOrWuISojd//eI1Nq15m3d7DmvIh7cPwMckEZsLzKKqqnvPMTVlZWeTk5JyxTps2bTCZTnb6mzVrFvfffz/5+flnPb7VasVqPbnoWGFhIbGxsRQUFBAQEHCuYQpx3uy2ch7/73haltjoU+6Ft6qjyCeGgXd/jF7W5BAerqAgled+vp4iWwVB5SYClFGkeo3GqTPw0hVdCfM3uztE0UwUFhZisVjO6fu7VulxWFgYYWH1tzaB2WzGbJZ/KMJ9DEYv7uvxMDlL3naV+ZceIWnTcjr1HeXGyIQ4u1/XvopVtVNhd5KuL+M4vxDkiKdD616ShAiPVW8/8VJTU9m6dSupqak4HA62bt3K1q1bKS6WuRmEZ4vpPgaHX5SmrHTDV6gO6SsiPNfxjK38kbkZu8OJ86+G7mC1LUWGNozuFOnm6ISoWb0lIs888ww9evRgxowZFBcX06NHD3r06MHGjRvr65RC1AlFp8en92RNmaH4GEc2/uKmiIQ4u/A9C7lRF4RXRWUvJwN6yo3X0SbMVyYwEx6t3hKRWbNmoapqlcewYcPq65RC1JkOvUdT5KPtLF2y4WuwSoue8EAZO1HSt9CxwszU/ACGlvkQpvaiQh/F2M7SGiI8m/S+E6Iaer0OU9+bXM9VVDbbMpm7/J9ujEqIaqgqbPkSgMxCKwYUEm2hFJquI8LiRY9YmcBMeDZJRISoQY9eA0jz68JRvY3Z/oUs8ClmwdFVHD++3d2hCXFSyu+Ql0xJhZ1iqx2Adb4jsenMjE2MRCcTmAkPJ4mIEDUwGXQ4e17Fd/7FZOgrP+ArHA6++f1lN0cmxF/sFbBtDlDZGgKQZwhll3cvLD5GBrQNcWd0QpwTSUSEOINRPXsTpnR1PVeBbfkHOHxQFsQT7udImg+l2ZTbHBSWVY7qWuM3FlXRM7pTJEaZ+0Y0AvIuFeIMfM0Gena+DzOVk/QFOQ20rxiEJaSPmyMTzV1JcQbPbfk/ljoLOVpUjgocNbUi2dQRb5OeoR3qb84nIeqSzPcrxFmM7daRLTtHYLGlkGWewn5DOMsPlXBFT+kEKNxnwdpXyXWW84Nahs6gMsjozTa/saAojIgPx9ukd3eIQpwTaRER4iyCfE0kdLqZQ773UWQIB2D53kzKKhxujkw0V1mZu/gtYz0AFXYnBYqDdT4hZBpjMOgVRiZEuDlCIc6dJCJCnIOxXaJQThl8UFbhYNW+LPcFJJq1XYcW4lCdqCpUOJzo0FFhvA6AQe3DsHgb3RyhEOdOEhEhzkGUxZseLbW3YhbvzsDmcLopItGcDev/EI8MfZVQXRiqCqFqIlZ9SxQFxiRKa4hoXCQREeIcnT5DZUGpjXWHcigpznBTRKI5i2pxEfn6GYQqV1FsmgJAn1bBhPt7uTkyIWpHEhEhzlHbMD86Rvq7nuucJSz6fQZPzbuKfQcWuDEy0RytOZBNcYWTQtNI7LrK1rpxnaPOspcQnkcSESFqYXyXyg96S8UiqHiKwxV/Uma38cPm93A67G6OTjQXDqfKwp3alrjEFhZahvi4KSIhzp8kIkLUQmJ0ALHBPihqIeVK5UyWVruTtPJs1m35yM3RieZiY0ou2cVWTdk4WdxONFKSiAhRC4qiMK5zJPmmSfiplbdpHE4VnCoVB5eDtcjNEYqm6OiRdez/8z1w2FBVlQWntYa0DvUl/pTbhkI0JpKICFFLvVsFExbgh1k3EYAYu5EbCy0Mc+hhy1dujk40NU6Hna/WvsTbe7/gyx8ms23rKtJySzV1xnWJQlFkcTvROEkiIkQt6XUKYxIjKTINpl9FF64p9sdcVjm3CIdWQMZOd4compDfN/6H1LJMANaVHeOjjY8RVrHWtT3C4kXPloFuik6ICyeJiBDn4aJ2ofh7Gdjvcz12xQzA8aLyyo1/flS5KqoQF6gg/zA/7//B9dzuVHE6vMg29nCVjescKa0holGTRESI82Ay6BjVKYJivYU1fmOAynlFym0OKD4OO793c4SiKfDf8QOX4If5r4/qCrsTL/2VqErlXCEWHyP924S4M0QhLpgkIkKcpxHx4XgZ9ezw7ku6sSUqkFn010iGPf+D3GS3xicauWNb0aWtY5jOn6f1kXR2emGxt6LIOMBVZXSnSIx6+RgXjZu8g4U4Tz4mAyMTwkFRWBYwCaeiJ6+kgnKbg93OEj5Z/ggOuUUjzoe9AjZ+6noaqBgYXhJGhel2V5mfl4FhHcPcEZ0QdUoSESEuwOjESMxGHXmGcDb4DKVUcfJheSbvObLYWnqUlevfdHeIojHa9QMUZ7qellbY+UU3gnKdxVU2NjESL6PeHdEJUackERHiAviZDYz6a8n1Tb5D+NbfySZKcagqAL8c+oXsrN3uDFE0NgVHKm/tnWKPPYrdXj1dz/28DAyPD2/oyISoF5KICHGBTrSKOBQDDsPVgEKFvXJVXpvq4PcN77g3QNFoqE4nZes/AOfJ5QJKbSpzGAunjIyR1hDRlEgiIsQFOrVVpMTYi1A1HpvdiUlVuC5qCJcNf8nNEYrGYt2Wj/hnxko2OktQ/2pV+904gFxDhKuOtIaIpkYSESHqwIlWEYAi003E2aPpFXw/Ay9+BcU70L3BiUahuCideXu+pggHs5w5/MeZRbLOh7nWvpp6Y6Q1RDQxkogIUQf8zAZGxlf+anXq/Dnm+zQrcluTUVDu5shEYzF39TOUOk+OskpSy5ln6INdMbnKfM0GRkhriGhiJBERoo6MToxwtYoAqCr8sv2YGyMSjYXDXoFJb9aUdfRtx67iPpqysZ2lNUQ0PZKICFFH/L2MrlaRE9YdypFWEXFWeoOJKWP/w+297sdHZ8JLZ0T1vllTR1pDRFMliYgQdehsrSLZWXv4ZeXTqE6nO8ITHq5r4rU8dcnnTOp0P9sy/TTbpDVENFWSiAhRh87UKrJx6yxeXnQHC1OXsmbju26KUHg6S2ArtuR10ZRJa4hoyiQREaKOVdcqMuuX+5m1/UPKnTYA5iZ9Q/qxje4KUXiwwzklbE3L15RJa4hoyiQREaKO+XsZGXFaq0hmSSucf80LAZUTnf1vzYtgtzZ0eMID2KzF2Aszqt3281ZtB2dpDRFNnSQiQtSDMae1ihSYRhJrONnc3lfx5cYKA2z+rzvCE242d8UTvPnz9Rzf+7OmXFpDRHMkiYgQ9aC6VpFDtqnEmUK5URfCjfoQvBQdHFgKh9e6KUrhDjt2fcvvmRtJdZbyyp8vs2bxI6jWYkBaQ0TzJImIEPXk9FYRJ95EBz9JX1OwtuKfH0FxVgNHJ9yhoCCVL7e+53pegcrc42vJL0yT1hDRbEkiIkQ9qa5VZEW6gZyEG7QVbaWw5t/gsCOaMFWl7M+P8Dtt6PZV8VMICkvgxy3SGiKaJ0lEhKhHp7eKqCrMyWoNcRdpK+YcgB3fNXB0okHtW0Rk1n4e00UwWKmcI6RrYDsG9LqLpIwith/J11SX1hDRXEgiIkQ98vcyulbmPWFzaj4HW00Gv5PlDlXlfzs/Z2/Sz6cfQjQF+Wmw5QsATIqOa/XB3OnThuuGvwqKwrcb0zTV/b2kNUQ0H5KICFHPxnaOxM/LoCn7ZlsO6sC/g85AjmrnbWcmi5wF/HfjvyksSKvhSKJRspVV3npzam+9db7oUfz8o9h4OI+U7BLNtsu6t5DWENFsSCIiRD3zMRmY2DVaU3Yws5jNJcFkJkzgZUcGh9TK+UQKHWV8sfxRmQK+qXA64Y93oOCItrzjeIjqht3hZO4m7bbwAC8Gtw9twCCFcC9JRIRoAMM6hhEeoF1d9ftNRwnqdA2t/VtqylOKUsk8vKohwxP15NCf/0E9ctoMuoEtodsUAFYmZZFVpJ3U7qpeLTDo5aNZNB/ybheiARj0Oq7oGaMpyyws5/eDuUwd+ToWgzcAwYqRh/o8TETr4e4IU9ShDVs+5c19XzPbmYf9xKy6Jj8Y/DAYTJRW2Pnfdu1ImbbhfvRsGeSGaIVwH0lEhGggveOCaBXqqyn7aetRjN6R3NjnIdoaAnhk+BtEJlzupghFXUlJWcVXOz8DYK1azDvOTIoABj8E/pWdlBfuzKC4XNtv5OpeMSiK0sDRCuFe9ZaIpKSkcMstt9C6dWu8vb1p27YtM2bMoKKior5OKYRHUxSFa3rHasqKyu0s3JlBx/YTuP+qH/GP6eOm6ERdsZYX8PEfz2FXT/bzOaha2dV+KER0AiCvpILFu45r9uvRMpD2Ef4NGqsQnqDeEpG9e/fidDr58MMP2bVrF2+99RYffPABTz75ZH2dUgiP1zHSn26xgZqyRbsyyC+tQDF5uycoUafMXhau6nwzRuXkqJchkf3o3+du1/Mftx7F5jiZqCiKwpW9tLfuhGguFFU9ZUnQevbaa6/x/vvvc+jQoXOqX1hYiMVioaCggICAgHqOToiGcSy/jGd+2smp//KGdgzjxgGtatwnI30LIcEdMJp9a6wjPEta2lo+/O1pwr2CuefSL9EbTAAcySvl2Z93aV7/YR3DmHqG11+IxqY239+GM26tYwUFBQQHB9e43Wq1YrWe7EFeWFjYEGEJ0aCiA70Z1C6U3/Znu8pW78tiVEIE0YFVW0V27ZnLpxv/TdegDkwb/zGKTrp2NQaxsQN4dMJM9HqjKwkBmLvpqCYJMRt1XNqthRsiFMIzNNgn2oEDB3jnnXe44447aqzz0ksvYbFYXI/Y2Nga6wrRmE3q0QKjXjv1++nzSQCsWvcmH2x4gwrVzsbc3cxfPaMhwxQXKMASi69fpOv5nvTCKlO5j0mMxOJjbODIhPActU5EHn/8cRRFOeNj7969mn2OHj3K2LFjufrqq7nttttqPPYTTzxBQUGB65GWJjNMiqYp0MfEmM7aqd+3puWz73iR63lBfgo/H5iHysmfzwtTl7Jz+xcNFqc4B07HOVVTVZXvT0s2A7yNjEmMrGEPIZqHWiciDz30EHv27Dnjo02bNq76x44dY/jw4QwcOJCPPvrojMc2m80EBARoHkI0VWMTo/A/ber37zamcaLbliWwFTf3fgCFk8M5eys+dNy1AI7vbtBYRfVysvfx7neXk7VvwVnrbkipOpX7pd2jZSp30ezVuo9IWFgYYWFh51T36NGjDB8+nF69ejFz5kx09XRv2+FwYLPZ6uXYouEZjUb0+qb/4ext0nNp92i+WpfqKjuUVcKmw3n0blXZlyox/gquLkjl26RvGKezMF4JQFEd8NsbcPFzYJGRFu6Sk72Pt5fcTa6thNfWvcAtufvp2O9vUM08ILZqpnKPsHgxuJ1M5S5EvXVWPXr0KMOGDSMuLo7XX3+drKws17bIyLppilRVlYyMDPLz8+vkeMJzBAYGEhkZ2eQndxrSPowlu4+TWXiyk/bczUfoHhvomuZ7SL/7iXNAXMrakztWFMPS52Dk05VThosGVZB3iP9bcg+5tsoWjlKc/Gff19zitNF94ENV6q9MyiK7WDuV+5U9Y2QqdyGox0RkyZIlHDhwgAMHDhATo/3VVlcjhk8kIeHh4fj4+DT5L63mQFVVSktLyczMBCAqKsrNEdUvg17HlT1jeH/lQVdZZqGVpXsyGdv5ZMIe1//vYLPB0VPWLbEWViYjI/4Bwa0bMuzmrSQHn9VvEGG3k3NKcZQ5mHZdrqtSvbDcxv+2aadybxfuR8+WgfUbpxCNRIPOI1JbZxqH7HA42LdvH+Hh4YSEhLgpQlFfcnJyyMzMpEOHDk3+No2qqvzr1z0cyjrZf8Bs1PHCpC4E+54c9omtHJY/DzkHNPvbDF5kDbiT6NgBDRVy81WaC0ufheLj2FSVT5zZ7FLLiDIH8fdxH+MfUHUY7qe/J/PHgWxN2RPj42kXLrOoiqarNvOINNp2wRN9Qnx8fNwciagPJ17X5tD3R1EUJvfV3l6x2px8/WeqtqLRq7L1I6yjq8imqnxsTeXNVY9xKHlZQ4TbfJXmwrJ/QnHl1OxGReFWXSjDfWJrTEL2HS+qkoT0aR0sSYgQp2i0icgJcjumaWpur2vbMD+GdNB2At98OK/KnBMYvWHYkxCRiF1V+dSZzW61nHKnjf/8/iz7DyxsuKCbkxNJSFG6pthoieHKCR9Vm4TYHU6+WHtYU+Zl1DO5j8yPJMSpGn0iIkRTcWWvGPxOG8771bpUrPbT5qkwesHQx1hmCWKnWuYqrlDtzFr7L2zF2sXUxIUpLzqOWk0Sgn9UZWdh78Bq91uy+zjH8ss0ZZf3aEGgj6na+kI0V5KICOEh/MyGKqvzZhdbmb89vWplg5kRF79JYsDJOXvM6Lg1cTpGv4iq9cV5KSo8yuvzp/NTwR5tJ3v/SBj5DHgHVbtfdrGVn0/roBob7MPw+PD6DFeIRkkSETeYPn16tTPSjh07tkHO/+yzz9K9e/cGOZeonYFtQ6osBb9wZ0aVX9YARrMvt034lG6B7TGhcE/CjbTueXNDhdrk5WTv4+1fbyWjIp+lzkJ+VPMrkxG/CBg5A3xqXjfr6/WpVNhPXV0XbhwQh17XvG45CnEuGnTRu/qkqipFVrtbY/A3G865b8PYsWOZOXOmpsxsNtdHWKIRURSFqQPiePbnXTidlb/AHU6VL9cd5pExHau8vwxGL24a9xHHDy6hRfxEd4TcJDkddt5b+neOVxS4ypY5i/DxCWLMWZKQLal5bE3L15QN7RhOmzC/+gpXiEatySQiRVY7D8zZ6tYY3prcnQCvc1u8ymw2Vzux28qVKxk9ejTLli1j8ODBALz66qu8/vrr7Nixg4iICBYuXMgLL7zAzp070ev1DBgwgLfffpu2bdu6jnPkyBEeeeQRFi1ahNVqJSEhgXfffZc9e/bw3HPPASc7hM6cOZPp06df4NWLutIi0JvRnSJYuDPDVZaUUcTaQzkMbFt1Jk6D0eusSYjDXqFZAVacmU5v4Nre9/HOH8+71voJMfrRd8S/wLfm6QLKbQ5mr9eOdvL3MnBlT1ldV4iayK0ZDzNs2DDuv/9+pk6dSkFBAVu2bOHpp5/mk08+ISKi8t5/SUkJDz74IBs3bmTZsmXodDouv/xynM7KpuDi4mKGDh3K0aNH+fnnn9m2bRuPPvooTqeTa6+9loceeojExETS09NJT0/n2muvdecli2pM7BatnUME+HZDGiXn0eq3fvNHvDHvanJzD5y9snDp0G4cV3W4CoAIk4W/X/wuQcFtz7jPL9vTyS2p0JRd0ycWH1OT+c0nRJ2Tfx1u8ssvv+Dnp22qffLJJ3nyySd54YUXWLJkCbfffjs7d+5k2rRpXHrppa56V155pWa/zz77jLCwMHbv3k3nzp2ZPXs2WVlZbNiwgeDgyibkdu3auer7+flhMBjqbKp9Ufe8jHqu69eS/yw/mTwUldv5YctRpvaPO+fjHDmyjjm7/otNdfDyr7dwc5+HiO946dl3FAAM6Xs/iqKnT5epeJ/hdgzA0fwyFu3K0JTFR/kzoI1MuCjEmUiLiJsMHz6crVu3ah533nknACaTia+++oq5c+dSXl7OW2+9pdl3//79TJkyhTZt2hAQEECrVq0ASE2tbBLeunUrPXr0cCUhonHq0TKI7rGBmrJVSZkcyio+p/1LS7P55LdnsKmVw39LnRW8u/5l0jd/Bp47oXKDK81LrvHvoeh0DOl331mTEFVV+WLtYVe/HgC9TuGG/nHNbk4cIWqrybSI+JsNvDW5u9tjOFe+vr6aVorT/fHHHwDk5uaSm5uLr6+va9vEiROJi4vj448/Jjo6GqfTSefOnamoqGwS9vb2Ps8rEJ5mSr+W7E4vdI3AUFX479rDPH1Jp7OOwCguyqhSNkjxJWrvIijKhgH3gMm3mj2bj3WbPmTu7i+4t8vtxHW/8byP88fBHPYfL9KUje0cSZRF/i0KcTZNpkVEURQCvIxufdTVL5+DBw/ywAMP8PHHH9OvXz+mTZvm6v+Rk5NDUlIS//jHPxg5ciQJCQnk5eVp9u/atStbt24lNze32uObTCYcDke124RnCfUzM7FbtKYsLbeUFXszz7pveERnHr1sNl0slQlvK8XEFbq/5r04ugkWPg55KXUdcqNgs5bw/ZIH+XLXLMpUBx/v+JjCtHXndaxiq51vN6ZpykL9zFzSNbqGPYQQp2oyiUhjY7VaycjI0Dyys7NxOBzccMMNjBkzhptuuomZM2eyfft23njjDQCCgoIICQnho48+4sCBAyxfvpwHH3xQc+wpU6YQGRnJpEmTWLNmDYcOHWLu3LmsXVu5jHyrVq1ITk5m69atZGdnY7Vaq8QnPMfoThFEBXppyuZtOUp+aUUNe5zk4xPK7RM/5/I2E7lFF4bx1GS5OBMW/wP2zgeHe4e+N6hjWzg4/15Wpq91FeWrdj757WnsRbWflXbupiMUl2v/ftf3b4nJIB+vQpwL+ZfiJgsXLiQqKkrzGDRoEC+++CKHDx/mww8/BCAqKoqPPvqIf/zjH2zbtg2dTsecOXPYtGkTnTt35oEHHuC1117THNtkMrF48WLCw8MZP348Xbp04eWXX3atYnvllVcyduxYhg8fTlhYGF9//XWDX784dwa9jqn9W2nKym0OZv2Rwrksnq3odIwc9CRBI2eA+bRVMB022Pxftv98B3uTfq7DqD1QSQ789gasfJn48lK6K9oFM2PDutQ4XXtNdh4tYPW+LE1Zz7ggusbU7jhCNGeKei6fZG5ypmWEy8vLSU5OpnXr1nh5edVwBNFYyetbVXXLyV/XryUjE2oxpXtJDvz+JuScHI1Tpjr5pyOdIhx0DWzHFQOeJDQsoa7Cdj+HHfYtgB3fgf1k61+eaud5RzoORWFyp6kM6HVnrQ5bWG7j2Z92UVB2coVos1HHC5O6VBl6LURzc6bv79NJi4gQjcTVvWMI8NZOmPftxjSOVjP9e418Q2DUc9B+tKtooVpAEZV9hrbnH+CFhbeSlbKyLkJ2u4IjGyr7wmz5UpOEAAQpBq4P7sajI96qdRKiqiqfr0nRJCEAV/SIkSREiFqSRESIRiLAy8jNF7XWlNkdKh+tOqhZ1+Ss9Abocwtc/E+yAqJY6dQOB+5otBDWon9dhOw2RYVH+eLXO3h2+X1k5R+qWsHkB/3upNeE92jRom+tj79yX1aVadw7t7AwMkEWtROitiQREaIR6RJjYVQn7a2YI3llzN18pPYHC+tIyLjXuSbhOnz1lesc6VG4svudYGy8t8NWr3+bf/48hfXZ27Gh8p0zT9uXpt0omPhvaDu8cjW6WjqWX8Y3f2pHyfh7Gbh5UGuZM0SI8yCJiBCNzJU9Y4gJ0s5PsXT3cXYcKahhj5rp9AYu6nMvMy7/nuFRAxkV2JHwjtWvW+OwV7D4txfIydl/XnE3FG+zP2XOk7dMdqvlbFPLIKgVjH4B+t4GZv+aD3AGNoeTj1YfwubQtkDdPKg1Fu9zW2dKCKEliYgQjYzJoOP2oW0x6rX/fD9bk0xhua2Gvc7MxyeUKy9+g4mXfFpjK8He/f/j5+T5PDt/Gu/Mu5ZN2z5HtZ99CHFD6544GT/9yRYds85AeYeLYcy/ILT9BR173uajpOWWaspGJITLKBkhLoAkIkI0Qi0Cvbm6d4ymrLDMxszfz21Ib410NX8krN1fObxXRSWpKJVl2z9D+fEu2DizwSZGKy3NZufu7/jkpxspPrK+2jpGow8DovoB0Cu4E89c8gX9+/4ddPoLOveuYwVV1pJpEeTN1b1iL+i4QjR3TWaKdyGamxHx4ew8Wsj2I/musu1H8lmRlMmI+FoM6T0HRYVH2Z6vXb13gOILFcWwb2HlwzsIAuMguDUHjEYsofGEhiagnCG5OSN7BeQfrhxqnHOQ74+uYmXZyb4wrXd8yciYftXuOqzX3QzsPJWw8MTzO/dpisptfPp7sqbMoFe4bXAbmbhMiAskiYgQjZSiKNw0qBUzftpF4SnDSL/dcISOkQG0CKy7dU6cThsXRfRiY9ZWypw2jCj0Uk5bp6Ysr/KRvpVP7UcpwoFZZ6CFdxhXd7mF2A4Tqhy3pDiDHUk/4VQdoKpEGP1oW1EBOQchPxXUk0sRBJ02uue37G2MKMlB8a26uq3F0rJuLpzKobqz1qRQUKq97XV1r1hig31q2EsIca4kEfEQw4YNo3v37vz73/8+r/2nT59Ofn4+P/74Y53GJTzbiSG9/166z1Vmczj5aNVBnprQqc5+rVsCW3HtmP/jcmsR23Z/Q96RP/Epqr5zbIHqcM1LYnXaOVSSjslZ/dpGObkH+HLXLNfzDooXf9dXPwS2paKdnyNbtbF352wS+v3tPK7o3K2SobpC1CtpU2yihg0bxv333+/uMEQDqGlI7w/nM6T3LExmf/r0uJXREz+CiW9D4uXgHaw9N9oOrCYUwiO6VXs8RdH220hVK2rs49ISEwqgV3S09onk2vgptO56w/lfzDk4ll/GHBmqK0S9khYRIZqAK3vGsDe9kCN5J2dZXbL7OInRFrrEWOrnpP6R0G0ydL0WCo9VdljNS6Es/U/88ooodpQDEKPzQrHEVHsI3WkdSMtxchw7kZwYCqtAQDSEtMMc0o7Hvf2IiOyBwVh3t51qUmGvfqjuTRfJUF0h6lLTTETKaz+fgovBGww1TNFcXghU82vNq3Yf9CUlJdx111388MMP+Pv78/DDD7u27d27l549e/LJJ59w3XXXAfDtt98ybdo0Nm3aRKdOnXA4HDzyyCN89tln6PV6brnlFs2vyOnTp7Nq1SpWrVrF22+/DUBycjKtWrWqVZyi8TAZdNw2pA3P/7Ibu+Pke+HD1Qd5YnxCnfYXqUJRwNKi8tHqInr3uJ5eTieFhakcSd8E1qIaR6yYjL609A5HARQUdIpCUVAikVG9IaQdBLUG08l+GC3q7yo0VFXl09+TqwzVHR4fTrfYwAaKQojmoWkmIj/cfv779r4ZOoypftv8Bys/VE933Te1OsUjjzzCqlWr+OmnnwgPD+fJJ59k8+bNdO/enfj4eF5//XXuvvtuBg0ahE6n48477+SVV16hU6dOALzxxhvMmjWLzz77jISEBN544w3mzZvHiBEjAHj77bfZt28fnTt35p///CcAYWFhtYpRND4xQT5c0zuW2etTXWVlFQ7+vWQfT01IINCn4dZAUXQ6LIGtsAS2OmO9sPBEHr36p4YJqha+33SEjSm5mrLoQG+u6S1DdYWoa00zEfFgxcXFfPrpp3z55ZeMHDkSgM8//5yYmJNN13fffTe//vorN9xwAyaTiT59+vC3v53skPfvf/+bJ554giuuuAKADz74gEWLFrm2WywWTCYTPj4+REZGNtCVCU8wIj6cpONFbErJc5XlllTw9rL9PDY2Hi/jhc2l0RysSMpk4U7tfCFeRj13DJWhukLUB0lEGtjBgwepqKigX7+T8x8EBwfTsWNHTb3PPvuMDh06oNPp2LVrl6tjXEFBAenp6Zr9DQYDvXv3vrCJrESToCgKtw5qQ17JXg5llbjKU3NK+XDVIe4d0Q69TjpZ1mT7kXy+WndYU6bTKdw9vC0xQTJUV4j6IOm9h9q2bRslJSWUlJSQnp7u7nBEI2Iy6PjbyPaEB5g15duP5DP7z1RJWGtwOKeED1Yd5PQ/z40D4kiMrqcOv0KIJtoicsVH57+v4Qyd+ia8SbWdVWuhbdu2GI1G1q9fT8uWlZMu5eXlsW/fPoYOHQpAbm4u06dP56mnniI9PZ3rr7+ezZs34+3tjcViISoqivXr1zNkyBAA7HY7mzZtomfPnq7zmEwmHI7q524QTV+Al5H7RnbgX7/uocRqd5Wv3JtJmJ+JsZ2j3Bid58kptvL20v1YbdoRMpd0i2Jwe+lfJUR9apqJSC1HsZz7cQMu+BB+fn7ccsstPPLII4SEhBAeHs5TTz2F7pRpsO+8805iY2P5xz/+gdVqpUePHjz88MO8++67ANx33328/PLLtG/fnvj4eN58803y8/M152nVqhXr168nJSUFPz8/goODNecQTV+kxYu/j2zHa4uSNCNpvtt4hBA/M31aBZ9h7+ajtMLOv5fup6BMO3PqgLYhTOreUON0hGi+5JvJDV577TUGDx7MxIkTGTVqFIMGDaJXr14A/Pe//+XXX3/liy++wGAw4Ovry5dffsnHH3/MggULAHjooYeYOnUq06ZNY8CAAfj7+3P55ZdrzvHwww+j1+vp1KkTYWFhpKamVolDNH3twv25dXCbKuWf/HaIA5nVjABrZuwOJ++uOMCx/DJNecdIf6YPbCWTlgnRABTVg28YFxYWYrFYKCgoICBA2xpRXl5OcnIyrVu3xsvLq4YjiMZKXt+6tWBHOt9v0s606ms28NSEBCICmuffV1VVPluTwh8HsjXlUYFePDEuAV9z02wwFqIhnOn7+3TSIiJEMzC2cyTDOmr7OpRY7fx76T4Ky2017NW0/bztWJUkJMDbyP2jOkgSIkQDkkREiGZAURSu6xdXZbr3zEIrby7eR35pRQ17Nj2qqvLrjnR+3npMU24y6Pj7yPaE+plr2FMIUR8kERGimdDrFO4c2rbK0vVpuaW8OH8PR0/rJ9EUOZwqX647zNzTblMpCtwxtC2tQ33dFJkQzZckIkI0I15GPfePak+Qr3a699ySCl76dQ97MwrdFFn9K7c5+M/yA6xMyqqy7bp+Lekua8gI4RaSiAjRzAT6mHh0TMcqE56VVTh4c/E+1h/KcVNk9aegzMZri5LYfiS/yrZr+sQyIj6i4YMSQgCSiAjRLIUHePHE+ATahGlvRTicKh+tPsSCHelNZgbWjIJy/jV/DynZJZpyg17hzmFtGZMo6zEJ4U6SiAjRTAV4GXl4TMdqb0l8v+kIX65Pxels3MnIgcwi/vXrHrKLrZpyH7OBh0d3lEndhPAA9ZqIXHrppbRs2RIvLy+ioqKYOnUqx44dO/uOQogGYTbouWd4O0YkhFfZtnJvJu+uOIDV3jiXCtiYkstri5I0U9wDhPiZeHJ8PO0j/N0UmRDiVPWaiAwfPpxvv/2WpKQk5s6dy8GDB7nqqqvq85RCiFrS6RSu69uSq3vHVtm2NS2f1xYmNbq5RhbvyuCDVQc1U9sDxIX48tT4TkRZzrCmlBCiQdVrIvLAAw/Qv39/4uLiGDhwII8//jjr1q3DZmtcH2oNYdiwYdx///0XdIzp06czadKkOonnTGbNmkVgYGC9n0c0HEVRGNs5kjuGtkWv005rnpxdwoyfdvHHgWyP7zeSWVjO20v3882GtCqr6HaJsfDo2I5YfIzuCU4IUa0Gmz4wNzeXr776ioEDB2I0ygdBfXj77bfr/IuiVatW3H///Zok6dprr2X8+PF1eh7hGfq2DsbibeSd5fspqzh5S6awzManvyezal8W1/eLo2WIzxmO0vCsdge/7khnwY4MHNX0axnaMYzr+8VVSbKEEO5X751VH3vsMXx9fQkJCSE1NZWffvqpxrpWq5XCwkLNQ5ydw+HA6XRisVgapKXC29ub8PCqfQpE09Ax0p8nxycQ4meqsu1AZjH//GUXX6w7TPFpfS/cQVVVNh3O5R/zdvLLtvRqk5AresYwtb8kIUJ4qlonIo8//jiKopzxsXfvXlf9Rx55hC1btrB48WL0ej033nhjjb/aX3rpJSwWi+sRG1v1nvW5KKoo0jwczpo7251e1+6s+cO1uKL4nOueSUlJCTfeeCN+fn5ERUXxxhtvaLZbrVYefvhhWrRoga+vL/369WPlypWu7Sdujfz888906tQJs9lMamqq5tbMRx99RHR0NE6nU3Psyy67jJtvvhmAgwcPctlllxEREYGfnx99+vRh6dKlrrrDhg3j8OHDPPDAA67X9tTzA+zbt6/Kaw7w1ltv0bZtW9fznTt3Mm7cOPz8/IiIiGDq1KlkZ2vX+RCeIzrQm6fGd6JrTGCVbapa2ZH1yR92sDIp020ja9ILynhzyT7eW3GQ3JKqU9QHeBu5a1hbJnSNklV0hfBgtb4189BDDzF9+vQz1mnT5uSy46GhoYSGhtKhQwcSEhKIjY1l3bp1DBgwoMp+TzzxBA8++KDreWFh4XklI0/89oTm+VP9niLKL6raujP+mEGF4+SH2CN9HiEuIK7aus+ve54S28m5CP7e8+90COpQ6/geeeQRVq1axU8//UR4eDhPPvkkmzdvpnv37gDce++97N69mzlz5hAdHc28efMYO3YsO3bsoH379gCUlpbyyiuv8MknnxASElKlheLqq6/mb3/7GytWrGDkyJFA5e2xhQsX8uuvvwJQXFzM+PHjefHFFzGbzfz3v/9l4sSJJCUl0bJlS3744Qe6devG7bffzm233VbttXTo0IHevXvz1Vdf8fzzz7vKv/rqK6677joA8vPzGTFiBLfeeitvvfUWZWVlPPbYY1xzzTUsX7681n8/0TAsPkbuG9WebWn5zNmQSmahdghsidXOF2sP89v+bK7r15K2YX4NEle5zcHP246xZPfxapMgRVEYlRDOpd2j8THJ4nVCeLpa/ysNCwsjLCzs7BWrceLXudVqrXa72WzGbG7aC04VFxfz6aef8uWXX7oShM8//5yYmBgAUlNTmTlzJqmpqURHRwPw8MMPs3DhQmbOnMm//vUvAGw2G++99x7dunWr9jxBQUGMGzeO2bNnu87z/fffExoayvDhwwHo1q2bZv/nn3+eefPm8fPPP3PvvfcSHByMXq/H39+fyMiaJ326/vrr+c9//uNKRPbt28emTZv48ssvAfjPf/5Djx49XLEDfPbZZ8TGxrJv3z46dKh9MicaTrfYQBKiAli8O4NftqVjc2hb2VKyS/jX/D10iw2kR8tAurYIrPMOoXaHk/2ZxWw/ks/6Q7kUlFXf4T0+yp/r+sXRIlBGxQjRWNTbz4X169ezYcMGBg0aRFBQEAcPHuTpp5+mbdu21baGNBcHDx6koqKCfv36ucqCg4Pp2LEjADt27MDhcFT5crZarYSEhLiem0wmunbtesZzXX/99dx222289957mM1mvvrqKyZPnoxOV3lHrri4mGeffZb58+eTnp6O3W6nrKyM1NTUWl3T5MmTefjhh1m3bh39+/fnq6++omfPnsTHxwOwbds2VqxYgZ9f1V/MBw8elESkETAZdFzSNZoBbUL4ZmMam1LyqtTZlpbPtrR8AFqH+tI1NpDuMYHEBnuf162RYqudHUcK2HYkn51HCzSdZ08X5Gvi2j6x9I4LktswQjQy9ZaI+Pj48MMPPzBjxgxKSkqIiopi7Nix/OMf/2jyrR4Xori4GL1ez6ZNm9Dr9Zptp36Re3uf/cN94sSJqKrK/Pnz6dOnD7/99htvvfWWa/vDDz/MkiVLeP3112nXrh3e3t5cddVVVFTUbkn4yMhIRowYwezZs+nfvz+zZ8/mrrvu0lzTxIkTeeWVV6rsGxVV/S0z4ZlC/MzcPawdu44VMHt9KhkF5dXWS84uITm7hJ+2HCXI10S3GAvdYgMJ9Tvzv32bw8me9EK2phVwILOoyhDc0+l1CmMSI5nQNQovo/7MlYUQHqneEpEuXbq47f7/S4Nf0jz3MdQ81PC5gc9pnnsbam7Sfbr/06io51S3Jm3btsVoNLJ+/XpatmwJQF5eHvv27WPo0KH06NEDh8NBZmYmgwcPrvXxT+Xl5cUVV1zBV199xYEDB+jYsSM9e/Z0bV+zZg3Tp0/n8ssvByoThpSUFM0xTCYTDsfZZ9a8/vrrefTRR5kyZQqHDh1i8uTJrm09e/Zk7ty5tGrVCoNB7tk3BYnRFp67NJGlezL537ZjlNtqfo/klVSwMimr2lVvL0TnFham9G1JpMWrTo8rhGhYTfJbwd907lM316aun+nCO+P5+flxyy238Mgjj7g6mT711FOu2yUdOnTg+uuv58Ybb+SNN96gR48eZGVlsWzZMrp27cqECRNqdb7rr7+eSy65hF27dnHDDTdotrVv354ffviBiRMnoigKTz/9dJVRNq1atWL16tVMnjwZs9lMaGhotee54ooruOuuu7jrrrsYPny4q38LwD333MPHH3/MlClTePTRRwkODubAgQPMmTOHTz75pErLj2gcDHodYztHMrh9KNvS8tl6JJ9dRwvPmJRcqLgQX7rFWugRG+Rxc5kIIc5Pk0xEPN1rr73mul3h7+/PQw89REFBgWv7zJkzeeGFF3jooYc4evQooaGh9O/fn0suuaTW5xoxYgTBwcEkJSW5RrGc8Oabb3LzzTczcOBAQkNDeeyxx6rM3fLPf/6TO+64g7Zt22K1Wmsceu3v78/EiRP59ttv+eyzzzTboqOjWbNmDY899hijR4/GarUSFxfH2LFjXQmYaLx8zQYGtgtlYLtQ7A4nSceL2H6kgG1p+WQVVd8x/VwZ9To6RQfQNcZCt5hAgnyrzm0ihGjcFNWD52wuLCzEYrFQUFBAQECAZlt5eTnJycm0bt0aLy9pmm1q5PVt/FRVJb2g3NVakpxVUu2EY6cL8jW5Eo/4KH/MBmkxE6KxOdP39+mkRUQIUS8URSE60JvoQG/GdYk65+UHZNSLEM2LJCJCiAYhCYYQojpyg14IIYQQbiOJiBBCCCHcRhIRIYQQQrhNo09EPHjQj7gA8roKIUTz0GgTEaOxclGt0tJSN0ci6sOJ1/XE6yyEEKJparSjZvR6PYGBgWRmZgKVa9tIr/zGT1VVSktLyczMJDAwUGZdFUKIJq7RJiKAa2n6E8mIaDoCAwNdr68QQoimq1EnIoqiEBUVRXh4ODabzd3hiDpiNBqlJUQIIZqJRp2InKDX6+WLSwghhGiEGm1nVSGEEEI0fpKICCGEEMJtJBERQgghhNt4dB+RE5NaFRYWujkSIYQQQpyrE9/b5zI5pUcnIkVFRQDExsa6ORIhhBBC1FZRUREWi+WMdRTVg+fSdjqdHDt2DH9//zqfrKywsJDY2FjS0tIICAio02N7Arm+xq+pX2NTvz5o+tco19f41dc1qqpKUVER0dHR6HRn7gXi0S0iOp2OmJiYej1HQEBAk32DgVxfU9DUr7GpXx80/WuU62v86uMaz9YScoJ0VhVCCCGE20giIoQQQgi3abaJiNlsZsaMGZjNZneHUi/k+hq/pn6NTf36oOlfo1xf4+cJ1+jRnVWFEEII0bQ12xYRIYQQQrifJCJCCCGEcBtJRIQQQgjhNpKICCGEEMJtmkwi8u6779KqVSu8vLzo168ff/755xnrf/fdd8THx+Pl5UWXLl349ddfNdtVVeWZZ54hKioKb29vRo0axf79++vzEs6qNtf48ccfM3jwYIKCgggKCmLUqFFV6k+fPh1FUTSPsWPH1vdl1Kg21zdr1qwqsXt5eWnqeNprWJvrGzZsWJXrUxSFCRMmuOp40uu3evVqJk6cSHR0NIqi8OOPP551n5UrV9KzZ0/MZjPt2rVj1qxZVerU9t91fartNf7www9cfPHFhIWFERAQwIABA1i0aJGmzrPPPlvlNYyPj6/Hq6hZba9v5cqV1b5HMzIyNPUa82tY3b8xRVFITEx01fGU1/Cll16iT58++Pv7Ex4ezqRJk0hKSjrrfp7wXdgkEpFvvvmGBx98kBkzZrB582a6devGmDFjyMzMrLb+H3/8wZQpU7jlllvYsmULkyZNYtKkSezcudNV59VXX+X//u//+OCDD1i/fj2+vr6MGTOG8vLyhrosjdpe48qVK5kyZQorVqxg7dq1xMbGMnr0aI4ePaqpN3bsWNLT012Pr7/+uiEup4raXh9UzgR4auyHDx/WbPek17C21/fDDz9orm3nzp3o9XquvvpqTT1Pef1KSkro1q0b77777jnVT05OZsKECQwfPpytW7dy//33c+utt2q+qM/nPVGfanuNq1ev5uKLL+bXX39l06ZNDB8+nIkTJ7JlyxZNvcTERM1r+Pvvv9dH+GdV2+s7ISkpSRN/eHi4a1tjfw3ffvttzbWlpaURHBxc5d+hJ7yGq1at4p577mHdunUsWbIEm83G6NGjKSkpqXEfj/kuVJuAvn37qvfcc4/rucPhUKOjo9WXXnqp2vrXXHONOmHCBE1Zv3791DvuuENVVVV1Op1qZGSk+tprr7m25+fnq2azWf3666/r4QrOrrbXeDq73a76+/urn3/+uats2rRp6mWXXVbXoZ6X2l7fzJkzVYvFUuPxPO01vNDX76233lL9/f3V4uJiV5knvX6nAtR58+adsc6jjz6qJiYmasquvfZadcyYMa7nF/o3q0/nco3V6dSpk/rcc8+5ns+YMUPt1q1b3QVWR87l+lasWKECal5eXo11mtprOG/ePFVRFDUlJcVV5qmvYWZmpgqoq1atqrGOp3wXNvoWkYqKCjZt2sSoUaNcZTqdjlGjRrF27dpq91m7dq2mPsCYMWNc9ZOTk8nIyNDUsVgs9OvXr8Zj1qfzucbTlZaWYrPZCA4O1pSvXLmS8PBwOnbsyF133UVOTk6dxn4uzvf6iouLiYuLIzY2lssuu4xdu3a5tnnSa1gXr9+nn37K5MmT8fX11ZR7wut3Ps72b7Au/maexul0UlRUVOXf4P79+4mOjqZNmzZcf/31pKamuinC89O9e3eioqK4+OKLWbNmjau8Kb6Gn376KaNGjSIuLk5T7omvYUFBAUCV99upPOW7sNEnItnZ2TgcDiIiIjTlERERVe5VnpCRkXHG+if+W5tj1qfzucbTPfbYY0RHR2veUGPHjuW///0vy5Yt45VXXmHVqlWMGzcOh8NRp/GfzflcX8eOHfnss8/46aef+PLLL3E6nQwcOJAjR44AnvUaXujr9+eff7Jz505uvfVWTbmnvH7no6Z/g4WFhZSVldXJe97TvP766xQXF3PNNde4yvr168esWbNYuHAh77//PsnJyQwePJiioiI3RnpuoqKi+OCDD5g7dy5z584lNjaWYcOGsXnzZqBuPrc8ybFjx1iwYEGVf4ee+Bo6nU7uv/9+LrroIjp37lxjPU/5LvTo1XdF3Xj55ZeZM2cOK1eu1HTonDx5suv/u3TpQteuXWnbti0rV65k5MiR7gj1nA0YMIABAwa4ng8cOJCEhAQ+/PBDnn/+eTdGVvc+/fRTunTpQt++fTXljfn1a25mz57Nc889x08//aTpQzFu3DjX/3ft2pV+/foRFxfHt99+yy233OKOUM9Zx44d6dixo+v5wIEDOXjwIG+99RZffPGFGyOrH59//jmBgYFMmjRJU+6Jr+E999zDzp073dbfqLYafYtIaGgoer2e48ePa8qPHz9OZGRktftERkaesf6J/9bmmPXpfK7xhNdff52XX36ZxYsX07Vr1zPWbdOmDaGhoRw4cOCCY66NC7m+E4xGIz169HDF7kmv4YVcX0lJCXPmzDmnDzR3vX7no6Z/gwEBAXh7e9fJe8JTzJkzh1tvvZVvv/22SjP46QIDA+nQoUOjeA2r07dvX1fsTek1VFWVzz77jKlTp2Iymc5Y192v4b333ssvv/zCihUriImJOWNdT/kubPSJiMlkolevXixbtsxV5nQ6WbZsmeYX86kGDBigqQ+wZMkSV/3WrVsTGRmpqVNYWMj69etrPGZ9Op9rhMrezs8//zwLFy6kd+/eZz3PkSNHyMnJISoqqk7iPlfne32ncjgc7NixwxW7J72GF3J93333HVarlRtuuOGs53HX63c+zvZvsC7eE57g66+/5qabbuLrr7/WDL2uSXFxMQcPHmwUr2F1tm7d6oq9qbyGUDki5cCBA+f0g8Bdr6Gqqtx7773MmzeP5cuX07p167Pu4zHfhXXW7dWN5syZo5rNZnXWrFnq7t271dtvv10NDAxUMzIyVFVV1alTp6qPP/64q/6aNWtUg8Ggvv766+qePXvUGTNmqEajUd2xY4erzssvv6wGBgaqP/30k7p9+3b1sssuU1u3bq2WlZU1+PWpau2v8eWXX1ZNJpP6/fffq+np6a5HUVGRqqqqWlRUpD788MPq2rVr1eTkZHXp0qVqz5491fbt26vl5eUef33PPfecumjRIvXgwYPqpk2b1MmTJ6teXl7qrl27XHU86TWs7fWdMGjQIPXaa6+tUu5pr19RUZG6ZcsWdcuWLSqgvvnmm+qWLVvUw4cPq6qqqo8//rg6depUV/1Dhw6pPj4+6iOPPKLu2bNHfffdd1W9Xq8uXLjQVedsf7OGVttr/Oqrr1SDwaC+++67mn+D+fn5rjoPPfSQunLlSjU5OVlds2aNOmrUKDU0NFTNzMz0+Ot766231B9//FHdv3+/umPHDvW+++5TdTqdunTpUledxv4annDDDTeo/fr1q/aYnvIa3nXXXarFYlFXrlypeb+Vlpa66njqd2GTSERUVVXfeecdtWXLlqrJZFL79u2rrlu3zrVt6NCh6rRp0zT1v/32W7VDhw6qyWRSExMT1fnz52u2O51O9emnn1YjIiJUs9msjhw5Uk1KSmqIS6lRba4xLi5OBao8ZsyYoaqqqpaWlqqjR49Ww8LCVKPRqMbFxam33Xab2z4gVLV213f//fe76kZERKjjx49XN2/erDmep72GtX2P7t27VwXUxYsXVzmWp71+J4Zynv44cU3Tpk1Thw4dWmWf7t27qyaTSW3Tpo068//buWPUxKIoAMNncHiojSAIERvBTtyDtq7A1iVYuAI7t2AT3Ii9exAUaxsXoHeKpAkJgTSeyHxfKXI518fDH3nX19dP6373mT3aT/c4Ho+/fX8pb0eWu91uqaqq9Hq9MpvNyuFweOzG3v10f+v1ugwGg1Kv10u73S6TyaTsdrtP6z7zNSzl7bhqo9Eom83myzV/yzX8al8R8eG++q3fhX/eNwAA8HBP/4wIAPC8hAgAkEaIAABphAgAkEaIAABphAgAkEaIAABphAgAkEaIACkmk0ksFovsMYBkQgQASOMv3oGHm8/nsd1uP7x2PB6j3+/nDASkESLAw12v15hOpzEajWK1WkVERKfTiVqtljwZ8Gh/swcA/j+tViuqqopmsxkvLy/Z4wCJPCMCAKQRIgBAGiECpKiqKm63W/YYQDIhAqTo9/ux3+/jdDrF5XKJ+/2ePRKQQIgAKZbLZdRqtRgOh9HpdOJ8PmePBCRwfBcASOMXEQAgjRABANIIEQAgjRABANIIEQAgjRABANIIEQAgjRABANIIEQAgjRABANIIEQAgjRABANL8A2rFz6MwwsqcAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kws = dict(kind='spline', s=0.7, periodic=True)\n",
    "x_dot_dxdt = dxdt(x, t, **kws)\n",
    "x_dot_derivative = pk.differentiation.Derivative(**kws)(x, t)\n",
    "\n",
    "plot_kws = dict(alpha=0.7, linewidth=3)\n",
    "plt.plot(t, x_dot, label='Exact', **plot_kws)\n",
    "plt.plot(t, x_dot_dxdt, '--', label='dxdt', **plot_kws)\n",
    "plt.plot(t, x_dot_derivative, ':', label='derivative', **plot_kws)\n",
    "plt.xlabel('t')\n",
    "plt.title('Equivalence of methods')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom differentiation method\n",
    "\n",
    "We also have the option of defining a fully custom differentiation function. Here we'll wrap numpy's `gradient` method. We can pass this method into  the `differentiator` argument of a `KoopmanContinuous` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-20T00:26:19.526100Z",
     "start_time": "2020-10-20T00:26:19.515760Z"
    }
   },
   "outputs": [],
   "source": [
    "from numpy import gradient\n",
    "\n",
    "def diff(x, t):\n",
    "    return gradient(x, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-20T00:26:20.349165Z",
     "start_time": "2020-10-20T00:26:20.336317Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>KoopmanContinuous(differentiator=&lt;function diff at 0x000001F2B61B2D40&gt;,\n",
       "                  observables=Identity(),\n",
       "                  regressor=PyDMDRegressor(regressor=&lt;pydmd.dmd.DMD object at 0x000001F29514A410&gt;))</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" ><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">KoopmanContinuous</label><div class=\"sk-toggleable__content\"><pre>KoopmanContinuous(differentiator=&lt;function diff at 0x000001F2B61B2D40&gt;,\n",
       "                  observables=Identity(),\n",
       "                  regressor=PyDMDRegressor(regressor=&lt;pydmd.dmd.DMD object at 0x000001F29514A410&gt;))</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" ><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">observables: Identity</label><div class=\"sk-toggleable__content\"><pre>Identity()</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-3\" type=\"checkbox\" ><label for=\"sk-estimator-id-3\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">Identity</label><div class=\"sk-toggleable__content\"><pre>Identity()</pre></div></div></div></div></div></div><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-4\" type=\"checkbox\" ><label for=\"sk-estimator-id-4\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">regressor: PyDMDRegressor</label><div class=\"sk-toggleable__content\"><pre>PyDMDRegressor(regressor=&lt;pydmd.dmd.DMD object at 0x000001F29514A410&gt;)</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-5\" type=\"checkbox\" ><label for=\"sk-estimator-id-5\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">PyDMDRegressor</label><div class=\"sk-toggleable__content\"><pre>PyDMDRegressor(regressor=&lt;pydmd.dmd.DMD object at 0x000001F29514A410&gt;)</pre></div></div></div></div></div></div></div></div></div></div>"
      ],
      "text/plain": [
       "KoopmanContinuous(differentiator=<function diff at 0x000001F2B61B2D40>,\n",
       "                  observables=Identity(),\n",
       "                  regressor=PyDMDRegressor(regressor=<pydmd.dmd.DMD object at 0x000001F29514A410>))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dmd = DMD(svd_rank=2)\n",
    "model = pk.KoopmanContinuous(differentiator=diff, regressor=dmd)\n",
    "model.fit(x, dt=t[1]-t[0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
